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ABSTRACT 

To facilitate geologic interpretation of satellite elevation 
potential field data, analysis techniques are developed and verified 
in the spherical domain that are commensurate with conventional flat 
earth methods of potential field interpretation. 

A powerful approach to the spherical earth problem relates potential 
field anomalies to a distribution of equivalent point sources by least 
squares matrix inversion. Linear transformations of the equivalent 
source field lead to corresponding geoidal anomalies, pseudo-anomalies, 
vector anomaly components, spatial derivatives, continuations, and dif- 
ferential magnetic pole reductions. 

Correlation of gravity and magnetic anomalies combined with other 
geological and geophysical data is useful for enhancing the quality 
and uniqueness of the geological interpretation of potential anomaly 
fields. Equivalent point source inversion is used to investigate 
qualitative visual spatial correlations of surface free-air gravity 
and POGO satellite magnetic anomalies and regional heat flow and tec- 
tonic data for North America and adjacent marine dreas. A more quan- 
titative analysis of the regional potential field anomaly correlations 
at satellite elevations also is considered utilizing Poisson's theorem 
in a moving-window linear regression analysis between suitable derivatives 
of the anomalous gravity and magnetic potentials. 

An inverse relationship is observed qualitatively and quantitatively 
between long-wavelength gravity and magnetic anomalies over continental 
terrane. Regional correlations of negative gravity and positive magnetic 
anomalies characterize areas of large relative crustal thickness and 
magnetization as exemplified by a prominent magnetic high which cor- 
responds to a gravity minima trend extending from the Anadarko Basin 
to the Cincinnati Arch. Negative magnetic and positive gravity anom- 
alies characterize thinner crust and regions of higher heat flow such 



as the Cordillera of North and Central America and the Yellowstone geo- 
thermal region* Over oceanic areas gravity and magnetic anomalies show 
limited correlation. 

The capacity to model potential field anomalies is essential to 
verifying and detailing the lithospheric interjfretation of these data. 
Gauss -Legendre quadrature integration is used to develop a versatile 
procedure for numerically modeling anomalous gravity and magnetic 
potentials and their spatial derivatives on a spherical earth due to 
a geologic body of arbitrary shape and physical properties. The 
procedure is well suited for satellite potential field modeling be- 
cause the accuracy and efficiency improves with increasing distance 
between source and observation points. 

The spherical earth modeling procedure is used to investigate 
Mississippi Embayment crustal structure synthesized from published 
surface wave dispersion , seismic refraction and Bouguer gravity 
anomaly studies of the region. Modeled gravity and magnetic anomaly 
comparisons with upward continued free-air gravity and reduced to 
the pole POGO satellite magnetic anomalies of the area support the 
failed-rift (aulacogen) hypothesis for the origin of the Mississippi 
Embayment. These results indicate that an aulacogen at satellite 
elevations is characterized by observable positive gravity and negative 
magnetic anomalies derived primarily from a high-density rift com- 
ponent of non-magnet ic lower crustal material. 

The procedures and results of this investigation provide a com- 
prehensive approach to the lithospheric analysis of potential field 
anomalies in the spherical domain and, hence, have widespread appli- 
cation in the analysis and design of satellite gravity and magnetic 
surveys for geological investigation. 
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I. INTRODUCTION 


A. Statement of the Problem 

The theory pf plate tectonics has revolutionized the geosciences 
by explaining paleo and contemporary geodynamics, as well as providing 
new insight into the processes which lead to mineral and petroleum 
deposits and unifying much of our geological knowledge into a coherent 
principle. The origin and much of the development of this concept has 
been based largely upon regional surveys and observations made in remote 
areas of the earth. The plate tectonics concept was triggered in fact 
by regional magnetic surveys of the ocean basins and earthquake epi- 
center and foci locations in areas generally remote from the scientific 
centers of the world. 

Today we have entered a period where technological developments 
offer the promise of another quantum jump in our knowledge of the earth. 
This technological development, satellite observations, already has had 
a tremendous impact upon surface-derived spectral information through 
the ERTS and LAND SAT programs. Now t„ese observations are being turned 
toward the study of the subsurface and in particular the earth's litho- 
sphere. 

Satellite observations have three principal advantages for the 
measurement of potential fields that are related to lithospheric varia- 
tions. First, data can be obtained from the entire earth - no region 
is too remote for observation - and patterns difficult to measure and 
perceive in surface or near-surface data are delineated by the broad 
data coverage. Second, satellite observations form a consistent data 
set free from non-uniformity caused by problems such as secular varia- 
tions in force fields and differences between continental and oceanic 
observations. Lastly, force fields observed from satellites, such as 


magnetic fields, are potential fields whose intensity is strongly 
dependent on the source to sensor distance. Thus, satellite-elevation 
observations are largely free from the near-surface geologic sources 
which tend to mask or distort the signatures of deep, broad litho- 
spheric variations. 

Of course, these advantages are negated in many applications 
because of the loss of resolving power with increased source to sensor 
distance. However, this is not necessarily the case for applications 
related to defining lithospheric provinces. On the contrary, the 
integration of the field expression of multiple sources permits the 
characterization of attributes of geologically significant litho- 
spheric regions measured in hundreds or even thousands of kilometers. 

These zones identified and characterized on a global basis should 
provide useful information for deciphering earth history including 
paleo and contemporary geodynamics, delineation of segments of the 
lithosphere into resource provinces, and for numerical modeling of 
lithospheric processes. 

The realization of the importance of these observations has led 
to increased availability of highly accurate data for analysis. How- 
ever, commensurate advances in procedures to analyze and interpret 
these data for lithologic information have not kept pace and are 
unavailable generally. The overall objective of this investigation is 
aimed toward rectifying this problem which currently limits the utility 
of satellite potential field data in lithospheric studies. 

B. Objectives 

Geologic interpretation of conventional low-level or ground-based 
gravity and magnetic survey data relies on a sophisticated arsenal of 
procedures to analyze the anomalies in the Cartesian plane. Conventional 
applications of procedures such as continuation, differentiation, magnetic 
reduction to the pole, etc., involve simplifying assumptions which gen- 
erally do not hold over large areas of a spherical earth that also exhibit 
considerable variations of geomagnetic field amplitude, declination and 
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inclination. Accordingly# attempts to directly adapt these procedures 
for analysis of regional-scale survey data registered in degrees of 
latitude and longitude are problematic and difficult to accomplish. 

An alternative approach is described and demonstrated in section 
II whereby many of the objectives of conventional processing technology 
can be achieved accurately and efficiently for gravity and magnetic 
anomalies registered in spherical coordinates. The procedure involves 
relating the potential anomaly ~ield to a spherical distribution of equi- 
valent point sources by least squares matrix inversion. The processing 
4 objective then can be computed in spherical coordinates directly from 

the equivalent point source distribution. 

However# the geologic interpretation of gravity and magnetic anomaly 
data# even wnen suitably processed# often leads to ambiguous results 
because of interfering anomalies# inadequate anomaly resolution and the 
multiplicity of equivalent anomaly sources. An important and time- 
honored technique for reducing this ambiguity is to analyze the cor- 
relations between the potential field anomalies along with other classes 
of geophysical data. Gravity and magnetic data are particularly suited 
for correlation analysis because anomalies in gravity and magnetic fields 
often are related genetically# although it is clear that the lack of 
anomaly correlation also can provide significant information for the 
geologic interpretation of these data. 

In section III# equivalent point source inversion is used to pre- 
pare long-wavelength gravity and magnetic anomaly fields over North 
America and adjacent oceanic regions for visual spatial correlation with 
| each other and regional heat flow and tectonic data. POGO satellite 

magnetic anomalies are differentially reduced to the pole so that the 
data uniformly reflect source magnetization characteristics. Surface free 
air gravity anomaly values also are upward continued to satellite ele- 
vations for comparison with the magnetic data and the long-wavelength# 
low order and degree spherical harmonic representation of heat flow. 

A more quantitative approach to the correlation of these potential 
field anomalies is investigated in section IV. Here# Poisson's theorem 
is used in a moving-window linear regression analysis performed between 


suitable derivatives of the anomalous gravity and magnetic potentials 
which are developed from equivalent point source inversion processing. 

The regression parameters generated by moving-window Poisson's analysis 
are useful for establishing patterns of correlation and non-correlation 
between the potential field anomalies. Por spatially correlating 
gravity and magnetic anomalies these results include an initial estimate 
of the ratio of physical property contrasts for the inferred source to 
facilitate further investigation and modeling. 

Modeling provides critical data for effective survey design and 
serves as an efficient approach for verifying and detailing anomaly 
interpretation. A technique is described and demonstrated in section V 
which utilizes Gauss-Legendre quadrature integration to numerically 
model gravity and magnetic anomalies in spherical coordinates due to a 
geologic source with arbitrary shape and physical properties. This 
procedure is used in section VI to model satellite-elevation gravity and 
magnetic signatures of the crustal structure of the Mississippi Ettbay- 
ment for comparison with the POGO satellite magnetic data and upward 
continued free-air gravity anomalies. 

In general, sections II through VI attempt to describe and verify 
techniques for generalized processing (e.g., differentiation, continuation, 
regridding, etc.), correlation analysis and modeling of potential field 
anomalies in the spherical domain that are commensurate with geophysical 
procedures of potential field interpretation enjoying widespread 
acceptance in conventional flat earth applications. These sections 
are documented so that they are individually complete and, hence, each 
section can be considered on 2m independent basis. However, because 
the spherical earth procedures are demonstrated principally by analyses 
of POGO satellite magnetic and surface free-air gravity anomaly data 
for North America, this documentation style includes some unavoidable 
repetition regarding the description of these data sets. 

Finally, in section VII the results of this investigation are 
reviewed in the context of utilizing satellite potential field mea- 
surements for lithospheric analysis. 
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II. SPHERICAL EARTH GRAVITY AMD MAGNETIC ANOMALY ANALYSIS 
BY EQUIVALENT POINT SOURCE INVERSION 


A. Abstract 

Artificial, earth-orbiting satellites are making increasingly 
available consistent, regional-scale potential field data which offer 
the premise of greatly enhancing our understanding of the structure, 
dynamics and geological history of the earth's crust and upper mantle. 

To realize this promise, procedures are required for lithospheric 
analysis of gravity and magnetic anomalies registered in spherical 
coordinates. An efficient and powerful approach to the spherical 
earth problem is to relate regional potential field anomalies to a 
spherical distribution of equivalent point sources by least squares 
matrix inversion. Simple linear transformations of the corresponding 
equivalent point source field yield a number of geophysically interesting 
quantities for the potential field anomalies because the anomalous 
potential and its spatial derivatives of all orders are linearly related. 
For gravity anomalies, for example, equivalent point source inversion 
can be used to compute corresponding geoidal deflections, anomaly 
continuations, vector anomaly components and derivatives at any point 
in source-free space. Similarly, for magnetic anomalies equivalent 
point source inversion can give least squares estimates at any point 
in source-free space of the vector anomaly components, anomaly contin- 
uations and spatial derivatives, as well as the magnetic anomalies 
differentially reduced to the pole. Pseudo-potential field anomalies 
also are available from equivalent point source inversion. A number 
of examples using l°-averaged surface free-air gravity anomalies and 
POGO satellite magnetometer data for the United States, Mexico and 
Central America illustrate the capabilities of the method. 
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B. Introduction 

RegionA 1-scale anomalies of the earth's magnetic and gravity 
fields provide important constraints for refining the basic model 
of plate tectonics and understanding intra-plate dynamics. These 
anomalies, especially when combined with other geophysical and 
geological data, may be anticipated to yield much information con- 
cerning regional physical property variations that are important for 
determining the structure, dynamics and geological history of the 
earth's crust and upper mantle. The realisation of the importance of 
these observations has led to increased availability of highly accurate 
regional-scale potential field data for analysis. Cosine nsur ate 
advances in data analysis for lithospheric interpretation, however, 
have not kept pace and generally are unavailable. 

Regional gravity and magnetic surveys represent spherically 
composited data for areas ranging in size from a few degrees to 
complete global coverage by satellite measurements. Typical processing 
objectives include data regridding, continuation, differentiation, 
and other wavelength filtering operations. For regional magnetic 
surveys it also is desirable from the standpoint of lithospheric 
interpretation to reduce the data for the variable geomagnetic field 
effects of inclination, declination and intensity. 

Spherical surface harmonic functions are a natural representa- 
tion for spherically registered gravity and magnetic anomaly fields. 

In general, spherical harmonics also are well suited for accomplishing 
the processing objectives described above because these operations 
involve only modifications of the geometric variables of the anomalies. 
However, the problem of representing potential field measurements over 
a spherical domain with an orthonormal set of spherical harmonic func- 
tions in general is a formidable task of confutation (e.g., Kaula, 1967), 
especially when it is desired to retain for analysis higher order and 
degree terms involving wavelengths of the order of a few hundred kilo- 
meters or less. 

When potential field anomalies are spherically registered at 
constant elevation, Fourier transform methods yield a simpler 
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numerical representation of the data which is suitable for strike* 
sensitive and wavelength filtering. These filters are applicable 
because the frequency response of the potential field anomaly is 
invariant for any uniform orthogonal translation of the degree* 
coordinates of the anomaly over the spherical domain. As an example , 
low-pass filtering a set of l°-averaged free-air gravity anomalies 
can be useful to help eliminate topographic mass contributions, where 
it is noted that topographic mass effects tend to be minimized when 

o 

free-air gravity anomalies are averaged over areas of the order 3 
and larger (Woollard, 1962; Strange and Woollard, 1964). 

Fourier transform methods applied in the spherical domain to 
discriminate anomaly source characteristics, however , normally neces- 
sitate regridding of the data to a Cartesian plane. The regridding 
procedure with all its attendant problems is required in such cases 
because the Fourier frequency characterization of an anomaly is 
variant for spherical translations of the anomaly source. Continua- 
tion and differentiation by Fourier transform methods also require 
Cartesian reregistration of the data to avoid distortions due to a 
Nyguist degree-wavelength which is variable in terms of the linear 
radial dimension of these operations. 

Another approach to spherical earth processing is to invert 
potential field anomalies on a spherical distribution of equivalent 
point sources. In contrast to processing by spherical harmonics 
or Fourier transform methods, equivalent point source inversion 
yields numerically efficient representations of potential field 
anomalies to any desired degree of wavelength resolution which can 
be directly processed in the spherical domain. 

In principle, equivalent point source inversion takes advantage 
of the ambiguity of source physical properties and distributions 
which is a feature inherent to all potential field anomalies. Every 
potential field observation, for example, represents the integrated 
effect of sources which may be characterized with equal facility by 
an infinite number of combinations of source physical properties and 
geometric parameters (e.g., Ramsey, 1940; Kellogg, 1954). The 
objective, then, is to relate a set of spherically registered potential 
field measurements to a distribution of equivalent sources with 
numerical properties which are amenable to spherical earth processing 
objectives. 
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Especially convenient fro* a numerical standpoint are equiva- 
lent point sources because the characteristics of their gravity and 
magnetic potentials are simplest to describe mathematically. As will 
be demo n strated, the point source potential and all its geophysically 
interesting derivatives can always be expressed as the product of an 
appropriate physical property value and a purely geometrical point 
source function which describes the inverse distance between Bource 
and observation points. By specifying the geometrical distribution 
of the equivalent point sources, the problem reduces to one of deter- 
mining a set of physical property values by matrix inversion techniques 
such that the potential field anosialies due to the equivalent point 
source distribution approximate the observed anomalies to a specified 
precision. Thus, once a suitable set of physical property values has 
bean found, the equivalent source field representation of the observed 
anomalies is calculated by sunning at each observation point the 
anomaly effect due to each of the equivalent point sources. To 
obtain an equivalent source field representation of the continuation, 
differentiation, etc. of the observed anomalies to the same degree of 
precision as reflected in the original approximation, the summation 
proc e dure is simply repeated using the appropriate geometrical point 
source function. 

To further illustrate the mechanism and versatility of equivalent 
point source inversion, it is desirable to develop explicitly the prop- 
erties of gravity point poles and magnetic point dipoles which facilitate 
analyses of potential field anomalies in the spherical domain for sub- 
surface information that are analogous to conventional flat earth methods 
of geophysical interpretation. 

C. The Gravity Point Pole Potential 

The gravity point pole is defined as a source of negligible dia- 
meter when compared to the distance between the source and observation 
points. Its scalar gravitational potential GP(R) at point with 
sp h erical coordinates (r ,0 ,$) located a distance R from the source point 
( r l'®l»^l* 9 iven by 
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GP(R) - 


G»Am 


where 


(XX-l) 


•8 3 2 

G ■ universal gravitational constant (>6.67*10 cm /(^*see ) ) , 

Am ■ mass contrast of the point pole, 

R ■ | r-r^ | • (r 2 +r 1 2 -2rr 1 cosfi ) ** , (II -2) 

r - radius vector directed from observation point (r,6,4) to 
earth's center, 

r. ■ radius vector directed from source point (r.,6.,9.) to earth's 
center, 

A » angle between r and ft at the observation point (r,6,+) , such that 
cos 6 - cosOcosd^+sinOsinOjCOsty-t^) , 

6,0, » co-latitude coordinates of observation and source points, 
respectively, and 

♦ ,6^ ■ longitude coordinates of observation and source points, respectively. 

The geometric characteristics of the point source potential are illus- 
trated in Figure II. 1. 

Because of the dependence of the potential on the inverse of 
the distance ft, the derivatives of the point pole potential GP(R) 
are linearly related by the relationship 
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This result provides a useful function for conside r in g derivatives of 
the potential. 

Now, the vector gravitational force field GF(R) at the distance 
R is given in spherical coordinates by the negative gradient of the 
scalar point pole potential as 


GF (ft) « -$(GP(R)) * (GP(R))ftR ■ ( -G |^{g)ftR)Am, where 


the gradient of R is 
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For o Mt of n- observation* (OBS) and o sot of p-point source geometric 
functions (Q) find o set of p -physical property values (AX) such that 

Qtr.ff.+ iri.Mi) * AX . 0BS(r,$,*) 


ure II. 1 Gravity and magnetic anomaly analysis in spherical coor- 
dinates by equivalent point source inversion (see text for 
details) . 
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< r> tg and d^ are the spherically orthogonal basis vectors . Accordingly* 
the gravitational force field at (r,6,9) is given by 

'■»<*» ■ !»«*» <k *t * I *> * I V’**- 

where the radial scalar component is 
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the horizontal 0-co- latitude scalar component is 

{-G |^(i) • IpffOd® - and <Xl-9) 

the horizontal 0- longitude scalar component is 

< fe 1 ? 1 * ‘sn, !?> »*- * { s jjf )4 “- ,II - 10) 

Also, because 

, 9CP(R>, 2 ,9GP(R>, 2 . , ,9r, 2 . ,1 9R, 2 . , l 9_R, 2 > 
‘“Ir -1 ■ (_ 5r ’ * ( IF* * ‘FI? 1 * We i? ’• 

the total horizontal component of the gravitational force field 
(I.t-7) can be expressed in terms of the radial component (II-8) as 
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Zn relating a set of gravity observations to equation (II-7), 
it is important to recognize the distinction between the total gravity 
anomaly and the gravity •f'-'rt due to local geologic mass variations. 

Xn geophysical practice, ite quantity which is considered is the anom- 
alous perturbation of the radial component of the earth's geocentric 
gravitational field. The use of this quantity, which commonly is called 
the gravity effect, Ag , for lithospheric analysis assumes the gravity 
anomaly due to geological mass variations is so small that perturbations 
of the earth's gravity field in an/ direction other than the normal 
are negligible. Thus, when relating an observed gravity anomaly value 
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to tho gravitational fore* field of » point pola, tha objective is to 
axprass tha obsarvad gravity affect as tha radial component of aquation 
CXX-7) which is givan by aquation (XX>8). 

For lithospheric interpretation * it frequently la desirable 
to analyse the apatlal derivatives of the gravity anomalies. Dif- 
ferentiation of a set of regional-scale gravity observations* for 
fcxuaple* can help identify the higher frequency crustal components 
of tha anomaly fiald. Tha derivatives also provide useful information 
for modeling the geometric characteristics of anomalous sources. 

To obtain equivalent source field representations of the spatial 
derivatives for a set of gravity observations it is necessary to 
develop the spatial derivatives of the gravity effect of the gravity 
point pole. The first order derivative components of the point pole 
gravity effect can be calculated from the negative gradient as 


■$<Ag) - -G {$£r))Ara - -G (A$(^r> ♦ (—jlhjAa 
R 3 R 3 R 


- {-G {A^ <iy) ♦ (ij) §|)fa)A», (11-12) 

R J R 

where the components of $R are given by equations (XX-4, S and 6) . 
Nov. to evaluate OA/dR) note that (<JA/3r) - OA/3R) (3R/3r) • 1 and, 
hence, 

OA/3R) - R/A. (11-13) 


Accordingly, equation (11-12) becomes 
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where the scalar first radial derivative of Ag is given by 
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the scalar first horizontal 9-co-latitude derivative of Ag is derived from 
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and the first horizontal 9-longitudinal derivative of Ag is obtained from 
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Also, by the sane argument which led to the development of equation 
(11-11) it can be shown that the total horizontal gradient of Ag 
may be expressed in terms of its radial derivative according to the 
relation 
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Another useful tool for lithospheric interpretation is the 
Laplacian which, in terms of the gravity effect of the point pole, 
is given by 
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Now, V (Ag) = 0 because Ag is a potential function and, hence, 
3r r 30 sin 0 3<J> 


where (3Ag/3r) is evaluated in (11-14) and the second radial derivative 
is given by 
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Equation (11-18) relates the horizontal gravity ancxnaly curvature 
to the radial curvature which is the Laplacian component normally 
computed for analysis. Because of the dependence of the Laplacian 
on first and second order derivatives it is especially sensitive to 
higher frequency anomalies. Hence, relating a set of regional scale 
gravity observations to the radial Laplacian component given by (11-18) 
can be useful for studying the shallower crustal anomaly components 
of the data. Also, the zero contour of equation (11-18) can be used 
to roughly delineate the lateral extent of the anomaly source because 
the edges of the source are approximated by the inflection points of 
the anomaly which correspond to zero curvature. From the sign con- 
vention used in equation (11-18) it is obvious that the zero contour 
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which delineates regions of positive or negative radial Laplacian 
values outlines , respectively, regions of positively or negatively 
contrasting densities. 

In consideration of the foregoing derivations it is clear that 
the results of this section are linearly related by geometry through 
equation (II-3) . Hence, once a set of gravity observations has been 
suitably related to a distribution of point poles according to equation 
(II-8) , then all the variables necessary for evaluating the gravita- 
tional potential, anomaly vector components and spatial derivatives, 
as well as the continuation of any of these quantities, also have 
been determined for the observation set. 

D. The Magnetic Point Dipole Potential 

The magnetic point dipole is defined as a dipole of negligible 
length when compared to the distance between the dipolar source and 
observation points. Its scalar magnetic potential MP(R) evaluated 
at the observation point (r ,0 ,<f>) located a distance R from the source 
point (r.,6-,4-) is given by 

MP (R) = Aj-V (1) = -A]-$(^), where (11-20) 

1 R R 

= the gradient operator in source point coordinates such that 

V- ^ = W observation point coordinates , 

1 R R 

R = given in (II-2) , and 

A j = the magnetic polarization contrast vector of the source point . 
The relationships between source point and observation point geo- 
metry are illustrated in Figure II. 1. 

As shown in Figure II. 2, the polarization contrast vector may be 
expressed in terms of its inclination, 1^, and declination, D^, and 
the spherically orthonormal basis vectors of the source point as 

A} = Aj*u^ * A j*{ (sinI 1 )^ r +(cosI 1 cosD 1 )e 0 +(cosI 1 sinD 1 )e^} , (11-21) 

where u^ * the unit polarization field vector. Here, it should be 
noted that if polarization is induced only, then the polarization con- 
trast is given by 
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A^ ■ Ak*?^, where (11-22) 

Ak ■ the magnetic susceptibility contrast, and 

* the geomagnetic field vector at the source point such that 

- p i*V 

According to equations (11-21) and (II-4, 5 and 6) the magnetic 
scalar potential of a point dipole reduces to 


MP(R) = -Aj-$(g) = ( )Aj 

= { -fR ( ^*“r^ )Aj 

“ ^2* {sinI i ( l7 ) + COSI i cosD i ( 7lf ) + COSI i sinD i ( 7ifee ff )} 

R 

= ^**{sinI^(A) + cosI^cosD^ (B) + cosI^sinD^ (C) } 

R 

= ^ *E. (11-23) 

R 3 

■f 4 

Now, the vector magnetic *.orce field, MF(R), at distance R due 
to the point dipole scalar potential is given by 


MF(R) = -V(MP(R)) = - { V(ic-)*E + ~ *7E>Aj 

R R 

- ‘ ( k ( 3 ,lE * E * 3 

R R 

= MF (r) e + MF(0)e o + MF(<(»)e. , where (11-24) 

r 0 

MF (r) MF (0) = { -?■ * J}A j , MF (^ ) = { ~ *j}Aj , and 

R R R 

, Rr cos6 - 

J = -T* sinI l +( — i t cosI i cosD 1 + 

Rr sin0. cos (<|i -«(>.) - 

(-i i ~ - f ) cos^sinD^ 

Equation (11-24) is useful for representing vector magnetic anomaly 
data. However, most magnetic data taken for geophysical applications 
represent scalar total anomaly intensity measurements which contain no 
information to discriminate the vector components. For geophysical 
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interpretation , then, the assumption is made that the magnetic anomaly 
due to geological magnetization variations is so small that perturbations 
of the geomagnetic field in any direction other than its direction are 
negligible. Kontis and Young (1964) have shown this assumption to be 
reasonably valid for anomalies with amplitudes less than about 10,000 
gamma* in the earth's magnetic field. Hence, for geophysical interpre- 
tation purposes, a total magnetic anomaly measurement is taken to repre- 
sent a perturbation of the geomagnetic field along its direction at the 
point of measurement. 

Thus, to represent a set of total intensity measurements in terms 
of the magnetic field of point dipoles it is necessary to compute equa- 
tion (11-24) in the direction of the geomagnetic field at the obser- 
vation points. This can be achieved by taking the vector dot product 

*► 

between equation (11-24) and the unit geomagnetic field vector, u, at 
the observation point* Accordingly, the total magnetic anomaly, AT, 
at R due to the point dipole is 

AT = u-MF(R) = {-$•*<£ ^ (|i)}Aj 

* MF (r) sinl + MF(0)cosIcosD + MF(<|>)cosIsinD, (11-25) 

where I and D are the inclination and declination, respectively , of 
the geomagnetic field at the observation point. In the application of 
equation (11-25) geomagnetic field models such as the IGRF-1965 (Cain 
et al,, 1967) commonly are used to obtain pertinent values of (I^,D^,F^) 
at the source point and (I,D) at the observation point. 

Consideration of equation (11-25) in its expanded form suggests 
that the derivation of derivatives for total anomaly fields is likely 
to yield fairly complicated algebraic expressions. Fortunately, these 
expressions can be considerably simplified by reducing the data to the 
pole which in effect is an adjustment of the anomalies to what they 
would be if radially (vertically) polarized by a magnetic field of 
uniform intensity. In conventional applications the spatial scale of 
the magnetic survey is small enough that a single value of geomagnetic 
field inclination, declination and amplitude may be assumed for the 
pole reduction. Over regional-scale areas, however, magnetic anomalies 
must be differentially reduced to the pole for values of geomagnetic 
inclination, declination and amplitude which are variable at all source 
and observation point locations. 

*) 1 gamma = 1 nanoTesla, nT 
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For total magnetic anomaly observations referenced to a distri- 
bution of equivalent point dipoles according to equation (11-25) , 
differential reduction to radial polarization is readily achieved by 
setting I«Ij-90°, D»D 1 =0° and using the normalized polarization A j ' 
given by Aj' ■ Aj (N/F^) where N is the desired polarizing amplitude. 
Accordingly* equation (11-25) reduced to the pole leads to a consider- 
ably simplified radial polarization anomaly AT(N) given by 

1 3A 2 

AT(N) * - - ~“g}Aj'. (11-26) 

R R 


An assumption which is implicit in the application of (IX-26) 
is* of course* that the magnetic anomalies contain no appreciable 
component of remanence. When this condition holds* equation (XI-26) 
can be used to adjust the magnetic observations so that they uniformly 
reflect anomaly source characteristics of geometry and magnetization. 
Because of this considerable practical significance to geophysical 
interpretation* the radially polarized anomaly* AT(N), in general is 
more suitable for derivative analysis than the total field anomaly* AT. 

The first order spatial derivatives of a set of vertically 
polarized magnetic anomalies may be obtained from the negative gradient 
of (11-26) which is 

2 2 

-$(AT(N)) - ^ {ij - ^}Aj' ■ $R {~ - )A j ' * 

R R R^ R 

where the scalar first radial derivative of AT(N) is given by 


3 AT (N) * 15A 3 9A , . . . 
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(11-27) 


the 


scalar first horizontal 6 -gradient component of AT (N) 


1 3AT(N) ,15A 2 B 9B,... 

R R 


is 

(11-28) 


and the scalar first horizontal ^-gradient component of AT(N) is 


1 SAT (N) r 15A 2 C 9C> . , 

’SIS"!* 


(11-29) 


Also, because AT(N) is a potential function its I^placian is 
zero* so that the radial Laplacian component of AT(N) may be equated 
with the horizontal component according to the relation 


(11-30) 
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where the second derivative of AT(N) is given by 
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E. Combined Gravity and Magnetic Point Source Potentials 

For a point mass which also is magnetized, the observation-to- 
source point geometry is invariant so that the magnetic and gravitational 
potentials of the point source can be compared directly. Consideration 
of equations (11-14 and 26), for example, shows that for a common point 
source the gravity effect and radially polarized magnetic anomaly are 
related by 


AT (N) 


G*Am 


-Aii 


G*Am 



(11-32) 


When the magnetic anomaly is calculated directly from the gravity 
potential as in equation (11-32) , the result commonly is called the 
pseudo-magnetic anomaly. The corresponding pseudo-gravity anomaly, 
on the other hand, is given by 


A9 


G*Am 

ty * 


AT(N) { 



(11-33) 


Pseudo-potential field calculations facilitate enhanced geologic 
interpretation of gravity and magnetic anomalies. The procedure , for 
example , is to relate an observed gravity anomaly to a set of point masses 
Am, by one of the least squares matrix inversion techniques described 
in the next section. Comparing an estimate for the bulk magnetization 
contrast-to-density contrast ratio of the anomalous body to the equiva~ 
lent point masses, Am, gives corresponding equivalent point magnetizations 
A j 1 , which can be used to calculate the pseudo-magnetic anomaly from 
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aquation (11-32) or, even more simply, from equation (11-26). A similar 
procedure can be used to develop from an observed magnetic anomaly the 
corresponding pseudo-gravity anomaly using equations (XX-33 or 8) . If 
the pseudo-potential field anomalies are similar to the observed 
potential field anomalies, common source structure can be inferred for 
both magnetic and gravity anomalies, thus, minimising the ambiguity 
inherent to the analysis of each of the anomalies by itself. Of course, 
dissimilarity between pseudo and observed anomalies also can lead to 
valid inferences regarding source lithology and structure. 

Equations (XX-32 and 33) are special consequences of a more 
general result noted by Poisson (1826) as a theorem which describes 
the relationship between the gravitational and magnetic potentials 
of any uniformly dense and magnetized body. In its geophysical 
context Poisson *8 theorem states that at an observation point the 
resultant magnetic anomaly is related to the first derivative of the 
gravity anomaly taken in the direction of magnetization by the 
magnetization contrast-to-density contrast ratio of the body. By 
equation (II-3) it is clear that Poisson's theorem also can be extended 
to relate the n-th magnetic anomaly derivative to the appropriate 
(n+l)-th gravity anomaly derivative. Equations (11-19 and 27), for 
example, combine to relate the second order radial gravity anomaly 
derivative to the radially polarized first order magnetic anomaly 
derivative for a common point source. In general, then, Poisson's 
theorem provides an important quantitative basis for comparing gravity 
and magnetic anomalies and, hence, for minimizing the ambiguity 
associated with the lithologic interpretation of these potential fields. 

F. The Inversion Problem 

As described previously, the measured potential field anomaly at 
any observation point is the summation at the observation point of the 
potential field effects due to each of the anomalous sources. Because 
of source ambiguity of potential fields it is clear that distributions 
of point sources exist for any measured anomaly such that the integrated 
anomalous effect of the point sources represents with equivalent facility 
the observed integrated effect of the true sources. The current 
objective, then, is directed toward resolving the problem of obtaining 
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for a given set of observed potential field anomalies, a distribution 

of equivalent point sources suitable for processing purposes. 

The effort begins from a consideration of equations (11*8 and 25) 

which relate equivalent point sources to measured anomalies of the 

earth's radial gravity field or total magnetic intensity field, 

respectively. Both of these equations may be generalized as shown 

in Figure II. 1, so that the incremental anomaly effect due to the p-th 

point source at the n*th observation point is given by a (R)*Ax , 

"P P 


where q^ (R) is the geometric point source function and Ax n is the 


appropriate physical property value. Accordingly, the relationship 


between a set of n*observations (obs. , obs. , . . . , obs ) and a set of 

l 2 n 


p-point sources can be expressed by the system of linear equations 

, n* A, i ♦ q n*‘*2 + — ♦ V 1 *, ■ ob *i 


q_,*Ax, ♦ q„*Ax_ + ... + q. *Ax 

^21 1 22 2 2p p 


obs. 


v ta i*v to j*-*v“p 


obs 


In matrix notation the linear system is more conveniently described by 


Q*AX = OBS, (11-34) 

where Q is an (r-row by p-col) matrix with elements given by the geo- 
metric point source functions (q^UO), AX is a (p by 1) column 


matrix containing the physical property values (Ax) , and OBS is an 


(n by 1) column matrix of the anomaly observations (obs ) . 

n 


Now, the geometric variables of the observation points are known, 
so that the complete determination of the elements of Q can be made by 
simply specifying the geometric distribution of the equivalent point 
sources. Hence, by pre-determining the equivalent point source loca- 
tions equation (11-34) can be solved for a set of physical property 
values AX by a number of well known matrix inversion techniques. 

For example, when p ■ n, Q is a square matrix and the exact 
solution given by 


AX * Q _1 *OBS 


(II-35) 
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can be obtained by well documented numerical methods (e.g., Carnahan, 
1969) if Q is not near-singular. Theoretically, Q has singular values 
only if the source point and observation point are specified at the 
same location, in which case the classical singularity of potential 
fields arises and the inversion effort is rendered meaningless. In 
practice, however, care must be taken to insure that the elements of Q 
are not made near-singular with respect to the working precision of the 
computer by improper scaling practice. An extreme example of bad 
scaling procedure is to specify source point locations so close to 
the observation point locations that the significant figures of the 
difference in their spatial location is in meters, whereas all other 
spatial variables of the Q elements have significant figures given in 
hundreds or thousands of kilometers. In such instances, the computer's 
ability to distinguish the spatial difference in source and observation 
point locations from its representation of zero may be so impaired 
that the resultant Q matrix can be considered singular within the 
computer's working precision. On the other hand, good scaling 
strategy as suggested by the example requires that all variables of any 
physical dimension of the problem be computed with roughly equivalent 
precision. With respect to the applications considered in this paper, 
experience has shown that proper scaling of the Q matrix in general 
is not difficult to achieve. 

2 2 

The solution given by (11-35) normally requires p = n machine 
storage locations for determining Q * . Hence, because of the large 
number of anomaly observations usually being considered, the exact 
solution often is impractical to compute even when the problem is 
suitably scaled. In most applications, then, computer storage con- 
straints require that problem (11-34) be set up so that the number (p) 
of equivalent point sources is substantially less than the number (n) 
of observations. 

When p<n, the solution to the overconstrained system of equation 
(11-34) is given by 

AX » (Q T Q)“ 1 *Q T *OBS. 


(11-36) 
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This solution is well known in the literature of nunerical methods 
as a least squares solution in the sense that it minimizes the 
Euclidean vector norm of the residual* |OBS - Q*Ax| . As s u mi ng 
the problem is properly scaled so that Q is not near-singular , 
procedures for obtaining (11-36) fall generally into two basic 
categories depending on how the values of the solution vector are 
to be used in subsequent analysis. 

For processing purposes such as continuation* differentiation* 
etc.* the solution AX merely represents a set of weighting coef- 
ficients which operate on suitably modified geometric point source 
functions to achieve a least squares approximation of the processed 
observations. The physical reality of the property values contained 
in AX are of limited consequence to these processing objectives and 
all that is required* simply, is a set of property values which 
provide a least squares fit between the observed anomalies and equi- 
valent point source anomalies. To obtain this type of solution* 
numerous methods are available in the general literature for application. 

The Choleski decomposition algorithm (Lawson and Hanson* 1974) * 

for example, can be recommended when it is desired to minimize machine 

computation time because the method yields a least squares solution 

set of weighting coefficients without the necessity of directly computing 

the product matrix Q^Q or its inverse. Choleski factorization is 

T 

applicable when Q is properly scaled because the product matrix Q Q 
is positive definite with Q as its Choleski factor. The application 
of Choleski decomposition requires at least (p ) storage allocations 
and* hence* the algorithm may not be optimal when machine storage is 
the primary constraint. 

Machine storage requirements can be reduced to roughly (p(p ♦ l)/2) 
allocations when the inverse (Q^Q) * is used to determine the solution. 
Because the product matrix Q T Q is symmetric* its inverse may be deter- 
mined using only those values contained in either the upper or lower 
triangular half of the product matrix. In general* then, the product 
matrix Q T Q can be operated on : n "packed" form where* for instance, the 
col'ums of the upper triangle of the product matrix are stored sequen- 
tially in a singly-dimensioned array of length (p(p + l)/2). Hence* 


i 
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to determine the solution (11-36) only about one half of the entries 

T *1 

of the product matrix need to be calculated to yield the inverse (Q Q) 
Although slower than the Choleski decomposition algorithm, the procedure 
minimises machine storage requirements for determining the least squares 
solution vector AX. 

The second major solution category occurs when i* is desired to 
constrain the values of the solution vector AX for subsequent analysis. 
This situation arises, for exanple, when physically realistic property 
values are required for modeling volume elements of the crust or some 
other geologic feature of the earth. Constrained solutions to equation 
(11-34) can be achieved by the eigenvalue-eigenvector matrix decom- 
position analysis suggested by Lanczos (1961) . 

Accordingly, the eigenvalue-eigenvector decomposition of Q is 
T 

given by Q • U*S*V to yield the solution vector 

AX » (V*8" 1 *U T )*OBS. (11-37) 

Here, U and V are orthogonal matrices with columns which are the eigen- 
vectors associated with the columns and rows of Q, respectively, and 
S is a diagonal matrix with elements which are the eigenvalues, X, of 

Q. The eigenvalues and eigenvectors are obtained from the square 

T -1 

symmetric product matrix Q Q and the inverse S is a diagonal matrix 
of eigenvalue inverses, 1/X. Jackson (1972) has shown that the Lanczos 
solution (XX-37) is a least squares solution for the ovarconstrained 
system (i.e., p<n) and, hence, is identical to the conventional solu- 
tion (XI-36) when Q is not near-singular. In contrast to conventional 
algorithms such as Choleski decomposition or packed form solutions, 
however, numerical procedures for obtaining the Lanczos solution require 

substantially greater machine computation time and storage of at least 
2 

(p + p(p + l)/2) allocations. 

The main advantage of the Lanczos solution is that eigenvalues 
contained in S can be manipulated to constrain the values of the solu- 
tion vector AX without altering its least squares property. Eliminating 
eigenvalues found near each end of the eigenvalue spectrum, for example, 
minimizes the occurrence of extreme property values in the least squares 
solution vector that may be difficult to rationalize physically. An 
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example of a more sophisticated procedure la presented by Vigneresse 
(1978) for tapering the eigenvalue spectrum to obtain a solution which 
Is constrained to lie within a restricted range of density contrasts 
In an application to gravity interpretation. 

The Lane ros procedure also provides a generalized least squares 
solution for the case where Q is near-singular and, hence. Where the 
conventional techniques are not applicable. Here, near-singularities 
of Q correspond to small eigenvalues which can be eliminated from the 
problem by comparing the eigenvalues to an acceptable threshold value. 

The surviving eigenvalues and associated eigenvectors are then com- 
bined as in (XX-37) to obtain a generalized solution which is a least 
r—4ares fit of the observed anomalies. Bowman et, al_. , (1979) present 
a fuller discussion of the procedure as applied to equivalent point 
source processing and some applications to gravity and magnetic 
anomaly inversion. 

Finally, the problem of specifying the geometric distribution 
of equivalent point sources for purposes of efficient inversion must 
be considered. Xn practice, machine computation time and storage 
requirements frequently dictate the source distribution characteristics 
for inversion. As discussed by Bowman at al. , (1979) , these requirments 
can be minimized by the strategic assignment of a point source to 
each feature of the observed data that is to be modeled. Here, the 
geometric characteristics of the strategic source distribution can be 
specified according to well known rules relating the depth and wave- 
length characteristics of point sources (e.g., Nettleton, 1942; Smellie, 

1956) . 

Another numerically convenient approach which is more automatic in 
its application, is to perform the inversion on a spherical grid of 
point sources at constant depth. For this application, the initial 
depth of the source grid should be on the order of the source grid 
spacing (e.g., Dampney, 1969; Ku, 1977). Xn practice, the source grid 
density of unknowns generally is specified according to available 
computer resources, so that adjustments of the equivalent source field 
to more precisely reflect wavelength characteristics of the observed \ 

data normally reduce to straightforward adjustments of the source grid 
depth. 
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Zt alto should bo no tod that edge offsets frequently oro oboorvod 
in radial continuations of tho equivalent sourco fluid when tho lutorul 
dinonsions of tho sourco grid aro roughly coincident with tbs o b s er v at ion 
grid. This result refleets the decreased energy density associated 
with point sources located near the edges of the distribution that 
tends to be emphasised by continuation. Edge effects are particularly 
troublescne when the observed data set is so large that processing has 
to be accomplished via subsets at one elevation which oust be rejoined 
at a new elevation. Experience has shown, however, that edge effects 
can be niniaised by including in the inversion a peripheral rind of 
ob s e rv ed data which extends laterally beyond the source grid a distance 
roughly equal to or greater than the radial distance between the source 
grid and the desired elevation of continuation. 

G. Examples of Geophysical Processing Applications 

To illustrate the equivalent point source inversion technique 

and its processing capabilities, regional gravity and magnetic anomaly 

data obtained from NASA-GSFC for the United States, Mexico and Central 

America are considered. Accordingly, Figure ZZ.3.A shows a stereo* 

o 

graphic equal -area polar (SEAP) projection of smoothed 1 -averaged 

o 

free-air gravity anomalies observed for the region .239.5-280.5) E 

longitude and (9.5-50.5)°N latitude. The observation grid consists of 

(42, 42) -values at 1° spacing where the relative elevation (Z) of the 

observations is 0.0 km. The amplitude range (AR) of the observations 

is between 41 mgal and -53 mgal and the amplitude mean (AM) is *0.61 aial. 

The least squares representation of the observations by »:p»ival* - 

point source inversion is given in Figure II. 3. B. Here, the observations 

are inverted on a source grid of (15,15)-values to obtain the least 

squares solution vector of density contrasts given by equation (XI-36) . 

o o 

The point source grid spans the region (240-280) E, (10-50) N at a depth 
of 400 km below the relative elevation of the observation grid. The 
point sources are uniformly spaced at 2.857° intervals over the ource 
grid. As shown by a consideration of the amplitude means of Fir^vs 
II. 3. A and 3.B, the m.tan error of the least squares fit of the 
observed gravity data is negligible. 







An obvious by-product of the inversion is a least squares density 

model for the observed gravity anomalies . By Bran's formula (Heiskanen 

and Moritz, 1967), if the anomalous potential of the model is known, 

then the geoid height anomaly due to the anomalous mass is given by 

2 

the ratio of anomalous potential to normal gravity (>9.8 m/sec ). A 
least squares estimate of the gravitational potential of the observed 
anomalies can be readily computed by integrating at each observation 
point the potentials due to each of the source points according to 
equation (II-l) . Hence, by dividing this result by the normal accelera- 
tion of gravity a least squares approximation of the geoidial anomalies 
can be obtained. Accordingly, Figure II. 4. A illustrates a least squares 
representation of the geoidal deflections due to the mass anomalies 
responsible for the observed free-air gravity anomaly data of Figure II. 3. A 
The geoidal deflections given in Figure II. 4. A are accurate on a relative 
basis. However, when compared to global-scale geoids, these anomalies 
must be adjusted by a nearly constant value that reflects the broad- 
scale geoidal deflection due to long-wavelength variations of the 
gravitational potential which are too large to be adequately sampled 
by the restricted spatial coverage considered in this illustration. 

Another processing example is given in Figure II. 4. B where the 
results of the inversion are used to obtain a least squares representa- 
tion of the observed gravity anomalies upward continued to a typical 
satellite elevation. In this application, the gravity effect of the 
equivalent point sources was recomputed at 450 km elevation using equa- 
tion (II-8) . Further processing applications are shown in Figures II.5.A 
and 5.B where the inversion results are used to obtain least squares 
representations at 450 km elevation of the radial first derivative and 
Laplacian component of the observed gravity data according to equations 
(11-14 and 18) , respectively. At 450 km elevation, it turns out that 
the radial Laplacian is nearly identical to the second radial derivative 
field of the gravity anomalies. Note also, that the anomalies given 
in Figures II. 4. A, 4.B, 5. A and 5.B all have now been reregistered to 
integral degree nodes of the spherical grid spanning the region (240- 
280) °E, ( 10-50) °N at l°-station intervals. 










In Figure IX. 6. A, an equivalent source dipole field representation 
of POGO satellite magnetometer observations conducted during 1968 is 
presented. The magnetic data were reduced in profile form according 
to the procedures described by Mayhew (1979) to yield 4271 anomaly 
values from 345 satellite tracks ranging in elevation between 240 km 
and 700 km over the study area. To reduce the total magnetic anomalies to 
a common elevation the profile data were inverted on a spherical grid 
of (15,15) spherical prisms* to obtain the least squares solution 
vector of dipole moments given in (11-36) . The equivalent prism source 
grid was located .t surface of the earth and uniformly spanned 
the region (240-280)°E, (10-50) °N at 2.857° intervals. The dipoles 
were oriented along the main geomagnetic field directions gi?en by 
the IGRF-1965 updated to 1968. The fit of the equivalent source field 
to the observed profile data is to 1 gamma standard deviation. 

Examples illustrating the general degree of this fit for an arbitrary [ j 

selection of profiles between 10°N and 50°N are given in Figure II. 7. J 

Recomputing the equivalent source field at 450 km elevation using i 

! 

equation (11-25) modified to include spherical prismatic volumes gives 
the least squares representation shown in Figure II. 6. A of variable 
elevation satellite magnetic anomaly data reduced to a common elevation 
of 450 km. S 

Assuming that the magnetization is predominantly by induction, the 
variable influence of the geomagnetic field on the total magnetic inten- ' 

sity anomalies recorded in Figure II. 6. A can be eliminated by differen- 
tially reducing the data to vertical polarization. This result is ; i 

achieved by recomputing the equivalent source field using equation 
(11-26) modified for spherical prisms. Accordingly, Figure II. 6. B 
illustrates the least squares representation of the total magnetic 
anomaly data of Figure II. 6. A differentially reduced to the pole using 
a normalized polarizing amplitude (N) of 60,000 gamma. 


*) To convert the spherical prismatic dipolar moment formulation of 
Mayhew (1979) to point dipole moments, it is necessary only to 
multiply the prismatic moments by the spherical volumes of the 
corresponding prisms. 
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A comparison of Figures II .6 .A and 6.B shows that anomalies 
in the southern portion of the radially polarized map generally are 
shifted differentially from 3° to as much as 5° to the north. Larger 
amplitudes also are observed for these anomalies because of the pos- 
itive difference between the scalar intensities of the radially polar- 
izing field and the IGRF-1965 updated to 1968 at low latitudes. 

Northern amplitudes, on the other hand, show only slight relative 
alterations which include principally northward displacements that 
rarely exceed 1° or 2°. The behavior of the radially polarized data, 
then, is consistent with that anticipated for the conventional pole 
reduction. 

In consideration of the foregoing examples it is clear that 
once a set of potential field anomalies have been related to a set 
of equivalent point sources by least squares matrix inversion, then 
a number of geophysically interesting characteristics of the observed 
data also are known with equivalent least squares precision. In 
general, these characteristics include the associated potentials, 
vector anomaly components and spatial derivatives, as well as con- 
tinuations of these elements. 

H. Conclusions 

Equivalent point source inversion is an efficient and versatile 
approach to spherical earth processing of regional potential field 
anomalies. The procedure is to relate observed anomalies to a set 
of equivalent point sources by a system of linear equations. By pre- 
determining the locations of the point sources, the problem simplifies 
to solving the linear equation system for a solution vector of physical 
property values. In practice, the system of equations normally is 
overconstrained in the sense that the number of point source physical 
property values of the solution is less than the number of observed 
anomaly values. From the general literature of numerical methods, 
however, several well-known matrix inversion methods are available 
to obtain the solution vector such that the resultant equivalent source 
field is a least squares representation of the observed data. 
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The anomalous potential and its spatial derivatives of all orders 
are linearly related so that when an equivalent point source least 
squares representation for any one of these elements has been determined, 
representations with equivalent least squares precision also are avail* 
able for the rest of the elements by simple linear transformations. 

This result has great practical significance because considerable geo- 
physical information is related to the anomalous potential and its 
spatial derivatives. Prom the least squares inversion of observed 
gravity anomalies, for example, geoidal anomalies, vector anomaly 
components and derivatives at any point in source-free spaco also may 
be computed with equivalent least squares precision. Similarly, the 
results of equivalent point source inversion of observed magnetic 
anomalies can be used to obtain least squares representations at any 
point in source-free space of the vector anomaly components, dif- 
ferential pole reductions and spatial derivatives. Pseudo-potential 
field calculations also can be readily performed using the equivalent 
point source inversion formulation. 

In consideration of these results, it is concluded that equivalent 
point source inversion represents a powerful and efficient method 
for spherical earth processing of regional-scale gravity and magnetic 
anomalies. Accordingly, the method has widespread application in the 
analysis and design of regional-scale gravity and magnetic surveys 
for lithospheric investigation. 
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III. QUALITATIVE CORRELATION ANALYSIS OF REGIONAL 
NORTH AMERICAN GRAVITY AND MAGNETIC ANOMALIES 
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A. Abstract 

Regional free-air gravity anomalies of North America and adjacent 
oceanic areas mapped by low-pass filtering and upward continuing to 4S0 
km elevation are compared to 450 km elevation POGO satellite magnetic 
anomalies differentially reduced to the pole, published 12th degree 
spherical harmonic representation of heat flow, and regional geologic 
information. In general, gravity anomalies of the continental areas 
parallel the prevailing north-south structural trends, while the pre- 
dominant oceanic anomalies strike east-west. The oceanic areas are 
more negative them the continental area except for the minimum over 
the Hudson Bay region. Free-air anomalies upward continued to 450 km 
elevation tend to correlate with the regional heat flow pattern. The 
positive free-air gravity anomaly associated with the Cordillera cor- 
responds to regional high heat flow values and the negative anomaly 
over Hudson Bay is correlative with a heat flow minimum. In addition, 
the generally negative anomalies of the oceanic areas correlate with 
low heat flow values. Comparison of gravity and satellite magnetic 
anomalies indicates little correlation over oceanic areas, although 
the oceanic areas have a more subdued anomaly pattern in both data sets. 
Over the continent, several prominent free-air gravity anomalies have 
an inverse correspondence to magnetic anomalies. Magnetic anomalies 
in contrast to gravity anomalies show limited correlation to the heat 
flow pattern. It is concluded that combined interpretation of the 
three geophysical data sets can decrease the ambiguity of the inter- 
pretation process because of the different parameters mapped and the 
variable range of the depth of the anomalous sources. 
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B. I ntroduction 

Long -wavelength geophysical anomalies covering several degrees 
of surface area have long been of interest in studying the mega- 
tectonic features of the earth. However, their potential usefulness 
has been hindered by limited areal coverage and the inherent ambi- 
guity in their interpretation. Recently, these problems have been 
alleviated somewhat by improved compositing and processing of data 
covering both continents and oceans and the availability of multi- 
data sets which When used in combined analysis tend to decrease 
possible alternative interpretations. Accordingly, regionally cor- 
relatable free-air gravity anomaly, differentially reduced to the 
pole satellite magnetic anomaly and regional heat flow maps of North 
and Central America and adjacent oceanic areas are presented and 
the correlation and potential significance of selected anomalies are 
discussed. 

Free-air gravity anomalies reflect differences in the mass per 
unit surface area of the subjacent earth which are calculated by 
comparing observed gravity with theoretical gravity values based on 
planetary (latitude) and elevation effects. Thus, free-air anomalies 
correspond to the combined gravity effect of topography and deeper 
geologic mass variations which may be difficult to separate for pur- 
poses of interpretation. This problem is compounded because major 
elevation variations occur as an effect of tectonism which is also 
the primary cause of geologic mass variations. Nonetheless, free- 
air anomalies sometimes are preferred to Bouguer anomalies for regional 
tectonic studies because the Bouguer topographic mass reduction 
eliminates the effect of large segments of mass that are involved in 
tectonism. 

When averaged over regional-scale areas, free-air anomalies are 
reduced in amplitude and tend toward a zero value, suggesting the 
mass per unit surface area down to some level in the earth is approxi- 
mately equal when considering large surface elements (Woollard, 1962) . 
This is strong direct support for the hypothesis of isostasy. Empir- 
ical evidence suggests that areas on the order of 3° and larger gen- 
erally tend to be characterized by free-air anomalies that approach 


zero (Woollen 6, 1962 t Strang* and Woo Hard, 1964), although tha sis* 
of a *p*cific area will ba a function of local strength character- 
lstlcs. Thus# regional scale variations of free-air gravity anomaly 
values on tha order of 3° or larger frequently are considered to 
reflect leostatic variations within the earth's crust and mantle. 
Averaged positive free-air anomalies corre s po n d to regions of mass 
excess or crustal loading# whereas negative anomalies correlate with 
regions of mass deficiency and# hence# the anomalies indicate the 
degree of isostatic equilibrum. 

A number of mechanisms have been postulated which may cause 
mass variations within the earth. These include density variations 
in the mantle due to convection# relatively static compositional or 
petrological variations, and the delayed response ot the earth to 
loading phenomena such as glaciation# de-glaciation# erosion and 
sediment deposition. Zn general, free-air gravity anomalies by 
themselves cannot be used to distinguish the causative mechanism 
due to the interpretations 1 ambiguity inherent to potential fields. 
However# when combined with other geophysical .and geological data# 
regional-scale free-air gravity anomalies may yield information con- 
cerning regional mass variations that are important to understanding 
the structure# dynamics and geological history of the earth's crust 
and upper mantle. Accordingly# regional anomalies in magnetics and 
heat flow also are presented and correlated with gravity anomalies 
to facilitate their interpretation. 

The correlation between these data sets is complicated because 
of the dissimilarity in causative characteristics and range of depth 
of origin. The relationship between rock properties which cause 
magnetic and gravity anomalies# for example # often is complex# 
non-linear and non-universal . Nonetheless# the correlation of these 
anomalies derived from continental crystalline crust is a standard 
approach to geological and geophysical interpretation (Hinze at al. , 
1975) . This is a consequence of the general correspondence between 
densities and magnetizations of specific crystalline lithologies 
(Nettleton and Elkins, 1944). However# there are many exceptions to 
the generalities because magnetic minerals are a minor constituent 
of rocks and are subject to alteration to non -magnetic forms and the 
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•f facts of strong remnant magnetisation. Similarly, hast flow will 
be reflected in temperatures that* in turn, alter the density of rocks 
by thermal expansion or contraction and control at least in part the 
effects of mineral alteration. Furthermore, heat flow and associated 
temperatures will control rock magnetisation through the Curie tempera- 
ture of ferrimagnetic minerals which is the temperature above which 
geologically significant magnetisations are eliminated. Thus, the three 
physical characteristics, density, magnetisation and temperature may 
be related within the crystalline crust of the continents, but they 
are subject to numerous complicating variations. 

The difference in source depth range further complicates the 
correlation of force fields, although this can be turned to the 
advantage of geologic interpretation. Gravity anomalies may originate 
from the surface to deep within the mantle or even the core, but the 
source of magnetic anomalies is limited by the depth of the Curie 
isotherm (“ 57S°C) which may be as deep as 70 km, but generally is 
at a depth of « 40 km. For source depth of heat flow. Pollack and 
Chapman (1977) have found that within continental regions 40% of the 
mean heat flux is derived from radiogenic sources within about the 
upper * 8 km of the crust. The remainder is obtained from radiogenic 
contributions of the lower crust and upper mantle and a deeper portion 
that enters the base of the lithosphere which varies in thickness 
from roughly 50 to 300 km. Thus, magnetic anomalies are limited to 
sources within the crust and uppermost mantle whilo gravity anomalies 
and variable surface heat flux have origins which may reach several 
hundred kilometers or more into the mantle. 

In general, then, correlation of these data sets, either direct 
or inverse, or even the lack of correlation can be useful in limiting 
the possible geologic interpretation of specific anomalies. Accord- 
ingly, the data sets provide different, but inter-related evidence 
which, when combined with surface geologic information, car. help to 
clarify the tectonic history of North America. 
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C. Data Sources and Preprocessing 

To investigate the regional wavelength features of North America'! 
free-air gravity anomaly field, a data set of l°-averaged surface 
values based on the GRS-1967 and IGSN-1971 was obtained from NASA's 
Goddard Space Flight Center for the region (214. 5-307. 5) °E longitude 
and (4. 5-58. 5) °N latitude. The observed data set consisted of 5170 
values (= 94 cols, 55 rows) at 1° intervals of which approximately 
600 were unassigned. Nearly all of the unassigned values were re- 
stricted to the regions (214. 5-242. 5) °E, (4.5-27.5)°N and (242.5- 
272.5)°E, (4. 5-14. 5) °N which are located in the Pacific Ocean of the 
southwest portion of the data set. In order to extend coverage to 
these regions, the missing values were interpolated from the surround- 
ing observations. The contoured stereographic equal-area polar 
(SEAP) projection in Figure III.1 illustrates the resultant data 
set with an amplitude range (AR) from 185 to -285 mgal and amplitude 
mean (AH) of -7.63 mgal. 

To minimize the effects of local topography as well as to enhance 
the longer wavelength characteristics of the data which are import- 
ant at satellite elevations, the gravity data were subjected to a 
low-pass filter that strongly attenuated anomalies of wavelengths 
less than 8°. The wavelength characteristics of the high-cut filter 
cure illustrated in Figure III. 2 and the filtered gravity anomalies 
are presented in Figure III. 3. Jlassically, unshaded anomalies 
(10 to -10 mgal) in Figure III. 3 correspond to regions of approxi- 
mate isostatic equilibrium. 

A further consideration in the design of the adopted filter was 

to obtain low-pass filtered anomalies which could be adequately re- 
o 

sampled at 2 intervals to optimize available computer resources for 
subsequent inversion processing. Accordingly, the smoothed gravity 
anomalies were resampled at 2° intervals to yield the set of (46,27) 
anomaly values for the region (216-306)°E, (6-58) °N illustrated in 
Figure III. 4. 

To effectively project the gravity anomalies to the satellite 
elevations of the magnetic anomalies, the resampled data were inverted 
on a spherical grid of (19,14) equivalent point gravity poles. The 
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equivalent source grid was located at a depth of 400 km below the 
relative elevation of the observation grid and spanned the region 
(216-306 )°e, (6-58) °N. The station spacing of the source grid was 
5° in longitude and 4° in latitude. The resultant equivalent source 
field approximation of the resampled data is given in Figure II I. 5. 

The field of the equivalent point source model was calculated at 
450 km elevation above sea level in Figure III. 6 to provide an esti- 
mate of the upward continued gravity anomalies for comparison with 
the satellite magnetic anomaly data. 

POGO satellite magnetometer observations obtained during 1968 
were processed using the procedures described by Mayhew (1979) to 
yield a magnetic anomaly map of North America at 450 km elevation. 
Briefly, data processing consisted of removal of the POGO 13th degree 
geomagnetic field from the original orbital profiled data. Next, 
data profiles were selected in the elevation range (240-700) km and a 
least squares quadratic function was fitted to the profiles and 
removed to account for DS variation arising from ring current flow 
in the magnetosphere during storm time. An inversion then was per- 
formed on a spherical surface grid of prismatic dipolar moments which 
were oriented in the local direction of the IGRF-1965 updated to 1968. 
The field of these dipolar moments was recomputed at 450 km elevation 
in the local direction of the IGRF-1965 updated to 1968 over the 
region (220-300) °E, ( 10-50) °N to obtain the total magnetic intensity 
anomaly map for North America given in Figure III. 7. 

Assuming that the anomalies are derived primarily by induction 
in the earth's present geomagnetic field, the magnetic map was dif- 
ferentially reduced to radial (vertical) polarization to remove dis- 
tortions incorporated in the data by variable inclination, declination 
and amplitude effects of the geomagnetic field. Differential reduction 
to the pole was achieved by recomputing the field of the equivalent 
sources assuming a normalizing induction field with an amplitude equal 
to 60,000 gammas and oriented radially at both source and observation 
grids. The differentially reduced to the pole magnetic anomaly map for 
North America is shown in Figure III. 8. 
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Finally, the regional heat flow nap for the study area (Chapaaa 
and Pollack, 1975) is shown in Figure XZZ.9 for comparison. This 
map is a smoothed, degree 12 spherical harmonic representation of 
observed heat flow which is supplemented in areas without observations 
by geologically predicted heat flow values. Major physiographic 
provinces and selected structures also are generalized for North 
America and adjacent marine areas in Figure ZZZ.9. 

D. Discussion 

Comparing Figures XXX. 3, 6, 8 and 9 shows a number of interesting 
associations which are outlined below. The possible geological signif- 
icance of these correlations also are considered. 

1. Free-Air Gravity Anomaly Field 

Although free-air anomalies near zero are generally considered 
indicative of isostatic equilibrium, it should be noted that care 
must be exercised with this interpretation because the averaged free- 
air gravity anomaly of the data set is approximately -7 mgal. This 
is reflected in the negative geoidal deflection over North America 
(Chapman and Talwani, 1979) and suggests a mass deficiency over the 
entire study area. Therefore, all observations regarding the isostatic 
state of selected portions of the study area should be considered on 
a relative basis. 

A study of the low-pass filtered surface free-air gravity anomalies 
(Figure XXI. 3) shows the continents characterized by near zero or local 
positive anomalies, whereas negative anomalies predominate in the 
oceanic regions. The major exception is the broad negative (-40 mgal) 
anomaly centered over Hudson Bay which has been interpreted to reflect 
incomplete isostatic adjustment of the crust following Pleistocene 
de-glaciation (Xnnes, et al. , 1968), although this explanation has 
been questioned (Jeffreys, 1975). 

Xn general, the Pacific Ocean is gravitationally more quiet 
than the Atlantic Ocean. The anomaly values in the Pacific Ocean 
increase to the west with increasing age of the oceanic crust, thus, 
suggesting a possible direct relationship between age and mass. 

However, oceanic ridge areas (e.g., the East Pacific and Juan deFuca 



Heat flow (Chapman and Pollack, (1975)), geological provinces and selected structures of 
North America. 




Ridges) of the Pacific are characterised by relative positive anomalies. 
In contrast, the western Atlantic along the eastern seaboard, despite 
its age and relative tectonic stability, is a broad negative gravity 
anomaly. Taiwan! and Le Pichon (1969) have postulated lateral inhomo- 
genieties in the upper mantle to explain this anomaly, although it is 
closely related to a deep oceanic plain. 

The largest variations of the entire free-air gravity data set 
are the negative (-60 mgal) anomalies associated with the Puerto Rico 
Trench. A common interpretation for these large negative free-air 
anomalies, which frequently are observed over oceanic trenches (e.g. , 
Vening Meinesz, 1964; Worzel, 1976), is that the earth's crust ie 
being subducted resulting in a mass deficiency. The Sigsbee Deep of 
the western portion of the Gulf of Mexico also is characterized by a 
large negative (-50 mgal) free-air gravity anomaly. This is an 
unusually large amplitude minimum for the observed bathymetry and, 
thus, raises questions regarding the nature of the crust beneath the 
Gulf of Mexico. 

Both oceans are characterized by pronounced east-west striking 
anomalies which contrast sharply to structurally compatible, general 
north-south orientations of the continental anomalies. Off the west 
coast of California and the Baja Peninsula, the east-west trending arms 
of the broad Pacific low (-20 mgal) show a particularly interesting 
correspondence to a number of major fracture zones. For instance, the 
northern-most west striking arm of this negative anomaly is closely 
bounded on the north by the Mendocino-Pioneer fracture zones and to the 
south by the Murray fracture zone. 

Within the continent, the Appalachians are associated with a 
weakly positive, free-air gravity anomaly which continues south across 
Florida and Cuba. The Mississippi Embayment also is affiliated with 
a weakly positive north-south trending linear anomaly. The mass excess 
of the Mississippi Embayment and Delta probably is due in part to an 
anomalously dense crust that may be a relic of the Mississippi Embayment 
aulacogen (Ervin and McGinnis, 1975) and not simply the added mass of 
the sediments. 

Within the eastern U.S., east of the Cordillera, the relation- 
ship between gravity anomalies and cratonic basins and uplifts is 
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inconsistent and poorly defined. The Ocala Uplift, for instance, 
is a relative positive anomaly while the Cincinnati Arch and Sioux 
Uplift are negative. The Anadarko Basin is associated with a negative 
anomaly in contrast to the Michigan, Illinois, Will is ton and Alberta 
Basins whicn are all characterised by relative positive free-air gravity 
anomalies. 

In general, positive gravity anomalies characterize the Cordillera 
of North and Central America. Maxima are located over the geothermal 
region of the Yellowstone area, the southern Rocky Mountains, the 
central plateau of Mexico, much of Central America south of Mexico 
and northern South America. Most of these areas are regional elevation 
highs, although the relationships are much too complex for the maxima 
to be explained simply by the excess topographic mass. A few examples 
illustrate this point with reference to Figures III. 10 and 11 which show, 
respectively, the l°-averaged topography of the study region and topo- 
graphic variations high-cut filtered for wavelengths roughly smaller 
than 8° using the filter of Figure III. 2. The Yellowstone and southern 
Rocky Mountain areas are topographic highs, but a topographic hign over 
the Idaho Batholith area which is larger than the Yellowstone area is 
not a maximum. The maximum over central Mexico is related to a high 
plateau area, buc is confined within higher mountainous belts. It is 
bordered on the north by the Cross Range Province which terminates the 
Basin and Range Province and on the south by the east-west Neovolcanic 
Province. In contrast, a maximum gravity anomaly occurs over the 
relatively low-lying Yucatan Peninsula. In southeast Mexico, the 
Cordillera gravity maxima are separated by a minimum which projects 
onto the northern extension of the Techuantepec Fracture Zone. 

The Coast Range of the U.S. and the Basin and Range Province of 
the U.S. and northern Mexico are generally free of gravity anomalies 
suggesting that they are approximately in isostatic equilibrium. 

However, the Basin and Range Province grades into a Volcanic Rift 
Province in northeastern California and southern Oregon and Idaho 
(Hamilton and Myers, 1966) which is a gravity maximum. This m a x i mum 
appears to be related to a mass excess which is related to Mesozoic 
crustal thinning and Cenozoic rifting and volcanism. This contrasts 
with the relative gravity minimum over the Columbia River Plateau 
basalts of Washington and adjacent states. 
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A study of Figure III. 6 shows that the major anomalies are still 
represented, although much of the anomaly detail has been lost in 
upward continuing the free-air gravity field to 450 km. The positive 
anomalies affiliated with the Cordilleran System of the U.S., southern 
Mexico and Central America have survived the upward continuation process, 
but with this exception the anomaly field is dominated by the negative 
anomalies centered over the Puerto Rico Trench, Hudson Bay, the Gulf 
of Mexico and the Pacific Plate. 

2 • Magnetic Anomaly Field 

The most obvious feature of the radially polarized magnetic 
anomaly map of Figure III. 8 is a prominent east-west high across 
the midcontinent which is breached by a negative anomaly in the 
Mississippi Embayment region (Mayhew, 1976; 1977). The closures of 
this magnetic high correspond to a trend of free-air gravity anomaly 
minima which extends from the Anadarko Basin of the eastern panhandle 
of Texas to the tri-state area of Indiana, Ohio and Kentucky (Figures 
III. 3 and 6). It is interesting to note that these areas of increased 
magnetic anomalies and decreased gravity anomalies also are shown by 
seismic evidence to be areas of thicker crust (Warren and Healy, 1973) 
as illustrated in Figure III. 12. Accordingly, the combination of 
negative gravity and positive magnetic anomalies can be interpreted 
to reflect displacement by anomalously thick crustal material of the 
denser and less magnetic material of the mantle. 

To the west, the large transcontinental magnetic high appears 
to overprint a northwest trend of negative magnetic anomalies which 
corresponds to the Cordillera and affiliated positive free-air 
gravity anomalies. Here the inverse correlation of gravity and 
magnetic anomalies may be related to increased heat flow from 
beneath the lithosphere which is likely to deflect the Curie isotherm 
to a shallower depth, thus, producing a region of relative net negative 
magnetization. Similarly, a negative magnetic closure over eastern 
Idaho which corresponds to a prominent positive free-air gravity anomaly, 
is interpreted as due to a lower net crustal magnetization related to 
a thinner magnetic crust over the path of the Yellowstone hot spot. 
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Other prominent features of Figure III. 8 include the magnetic 
high over the Anti 1 les-Bahamas Platform which extends to southern 
Florida and the broad negative magnetic anomaly that projects along 
the Gulf Coast of the U.S. from Louisiana westward across Mexico. The 
eastern component of the negative anomaly appears to correspond with 
the southern terminus of the weakly positive free-air gravity anomaly 
of the Mississippi Embayment, whereas the western component roughly 
correlates with a large positive gravity feature of the Cordillera. 

The magnetic anomalies of the western third of the continent 
(i.e. , the Cordillera) are relatively subdued in contrast to the char- 
acteristics of the eastern two-thirds of the continent. The western 
portion has relatively broad strong positive gravity anomalies and a 
regional heat flow maximum (Figure III. 9) with a shallower Curie iso- 
therm. Thus, the subdued magnetic anomaly pattern appears to be the 
result of a relatively thin and laterally homogeneous magnetic crust. 

In general, the correlation of magnetic and gravity anomalies 
over continental regions may be restricted because the gravity data con- 
tain long-wavelength effects of deep crustal or upper mantle sources 
which are too deep to be magnetic anomaly sources as a result of regimes 
where the temperature exceeds the Curie point of magnetite. The con- 
trasting north-south gravity trends and east-west magnetic trends 
over North America may be evidence of this effect. The correlation 
of anomalies also may be hindered by remanent magnetization effects 
which are at a large angle to the present field or may even be reversed. 
Finally, correlations may be limited by anomaly sources with incon- 
sistent density and magnetization distributions. However, even this 
information from restricted gravity and magnetic correlations can be 
useful as a guide to delineating and modeling sources by pointing out, 
for inst# je, remanent magnetization effects and possible depths of 
anomaly sources. 

3. Heat Flow Field 

Other than the previously noted correlation of the high heat 
flow region of western U.S. (Figure III. 9) with a subdued magnetic 
anomaly pattern (Figure III. 8), there is no general correlation of 
the two maps. However, in a regional sense, the direct correlation 
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of the heat flow map with the upward continued gravity anomaly map 

(Figure III. 6) is rather striking. This supports the suggestion 

that regional gravity and heat flow anomalies may have a deeper source 

than magnetic anomalies and, thus, may be more prone to correlation. 

The positive gravity anomalies of the Cordillera correlate with 

-2 

anomalously high heat flow values (>60 mWm ) , but the axis of the 
heat flow high lies west of the axis of the gravity anomaly. This is 
believed to be the result of the very broad wavelength of the heat 
flow high (12th degree spherical harmonic) which has superimposed the 
effect of the heat flow high of the magmatic region of western North 
America with the heat flow maxima associated with the Gorda and Juan 
deFuca Oceanic Ridges off the northwest coast of the U.S. In general, 
the negative gravity anomalies of the oceans correlate with low heat 
flow values, although an exception occurs for the observed heat flow 
maximum centered over the Baja Peninsula (Chapman and Pollack, 1975) 
which correlates with negative gravity anomaly values. The negative 
gravity anomalies of Hudson Bay and the western Atlantic Ocean also 
are correlative with relative heat flow minima. 

E. Conclusions 

Regional continental free-air gravity anomalies of North America 
parallel prevailing north-south structural trends, while the anomalies 
of the adjacent oceanic areas are more negative and subdued and strike 
roughly east-west. The mean free-air anomaly of the area is approxi- 
mately -7 mgal suggesting a broad general mass deficiency which is 
reflected in negative geoidal deflections. Western North America ana 
much of Central America are characterized by broad positive free-air 
gravity anomalies in contrast to more areally-limited maxima and minima 
of the eastern protion of the continent. Hudson Bay and western Gulf 
of Mexico are marked negative anomalies. Some regional gravity 
anomalies are related to surface elevation or bathymetry, but the 
number of exceptions are sufficient to suggest that the sources of 
the anomalies are largely deeper geologic mass variations. This is sub- 
stantiated by the correlation of several regional gravity and magnetic 
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anomalies, although there are exceptions wnich may have different 
explanations. However, it appears that several gravity anomalies 
are derived from beneath the magnetic crust. 

The magnetic anomaly map shows a more complex pattern with 
positive anomalies over the continents and Antilles-Bahamas Platform 
as compared to generally negative anomalies for the oceans. The 
entire area has an unexplained east-west striking bias to the magnetic 
anomalies. The general relationship between gravity and magnetic 
anomalies over the continents is inverse. Regions of large relative 
crustal thickness in particular appear to be characterized by negative 
free-air gravity and positive magnetic anomalies. Negative magnetic 
and positive gravity anomalies, on the other hand, may be related to 
regions of increased heat flow. 

Free-air gravity anomalies upward continued to 450 km directly 
correlate with a 12th degree spherical harmonic representation of 
the heat flow, but little correlation is observed between specific 
magnetic anomalies and this heat flow pattern. Nevertheless, high 
heat flow areas tend to have a subdued magnetic pattern as a result 
of the lower net crustal magnetization associated with shallow Curie 
isotherm regions. 

These broad conclusions suggest that combined interpretation of 
long-wavelength gravity, magnetic, and heat flow anomalies can 
significantly decrease the ambiguity of the interpretation process 
because these anomalies map not only different parameters, but their 
potential sources tend to cover different depth ranges within the earth. 
In general, gravity anomalies and heat flow are derived from both 
crustal and sub-crustal sources, while magnetic anomalies have prin- 
cipally crustal sources. Therefore, the mapping of long-wavelength 
magnetic anomalies provides a valuable complementary set of data to 
gravity and heat flow anomalies for the investigation of megatectonic 
features of the continents and oceans. 
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IV. QUANTITATIVE CORRELATION ANALYSIS OP REGIONAL 
NORTH AMERICAN GRAVITY AND MAGNETIC ANOMALIES 


A. Abstract 

To quantitatively investigate the regional correlations of 
gravity and magnetic anomalies over North America, Poisson's theorem 
is used in a moving-window linear regression analysis performed 
between suitable derivatives of the anomalous potential fields at 
450 km elevation* The results of moving-window Poisson's analysis 
are useful for establishing patterns of anomaly correlation and pro- 
vide initial physical property constraints to facilitate further 
investigation of correlating anomalies. The statistical correlation 
of regional gravity and magnetic anomalies over continental terrane 
is predominently inverse. In particular, regions of large relative 
crustal thickness are associated with negative gravity and positive 

magnetic anomalies, as well as exceptionally large magnetization^* 

-2 3 

(- 1.8x10 emu/cm ) . This characterization is exemplified by a 
prominent east-west magnetic high across the midcontinent which 
corresponds to a trend of gravity maxima that extends from western 
Texas to the Cincinnati Arch. Negative magnetic and positive gravity 
anomalies, on the other hand, are related to thinner crust and regions 
of higher heat flow. Typical examples include the Cordillera of North 
and Central America and the Yellowstone geothermal region. Over 
oceanic areas regional anomaly correlations are poorly developed 
probably because of decreased crustal magnetization. Statistically, 
however, these regions are characterized primarily by direct gravity 
and magnetic anomaly correspondences. 

*) 1 emu/cm 3 ■ 1 oersted ■ (10 3 /4ir)A/m 
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B. Introduction 

The utility of long -wavelength gravity anomalies covering several 
degrees of surface area, for geologic analysis has been recognized for 
many years. More recently magnetic anomalies which have significant 
energy in wavelengths of several hundred kilometers or more have been 
identified as originating from within the earth's crust and perhaps 
the uppermost mantle (Alldredge et al . , 1963) . The geological inter- 
pretation of these regional anomalies is hindered by the effect of 
anomaly superposition and the source ambiguity which is inherent to 
the analysis of potential fields. A common approach to minimizing 
these interpretations 1 limitations, especially in continental terrene, 
is to evaluate the correlation between anomalous gravity and magnetic 
potentials. 

The basis of correlation analysis is the hypothesis, commonly 
validated in continental areas, that variations in lithology and physical 
properties of the crystalline crust are reflected in both their bulk 
mineralogy, which controls rock density and, thus, gravity anomalies, 
and their minor mineral content, specifically magnetite, which is the 
primary cause of magnetic anomalies. Consideration of the general mag- 
netic and density properties of the crust and upper mantle and the 
geologic causes of their variations suggests that correlation between 
long-wavelength magnetic and gravity anomalies, both direct and inverse, 
are anticipated, but not on a universal basis. Seismic and other geo- 
physical and geologic evidence shows that the lithosphere is marked by 
vertical variations, principally by the Moho which separates the crust 
and mantle and by a much more poorly defined and less consistent boundary 
between the upper and lower crust. Increasing evidence points to lateral 
variations in the petrophysical and chemical attributes of these so 
called "layers." The vertical and lateral variations plus temperature 
perturbations are all potential sources of gravity and magnetic anomalies. 

Geophysical and geologic evidence clearly indicates that the den- 
sity of the crust increases from a mean value of approximately 2.7 gm/cm^ 

3 

at the surface to roughly the 2.9 to 3.0 gm/cm range at the base of 
the crust where there is a sharp increase to roughly 3.3 gm/cm^ in the 
upper mantle. The magnetic properties of these layers are much less 
precisely known, but generally the lower crust is believed to hava a 
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magnetization of roughly 5x10 emu/cm made up equally of induced 
and viscous remanent magnetization, which is several orders of mag- 
nitude greater than the mean magnetization of the upper crust. 

Wasilewski at al . , (1979) present a strong case for a non-ferrimagnetic 
upper mantle. As a result for a constant thickness crust, we 
anticipate a direct correlation between gravity and magnetic anomalies 
as the thickness of the layers of the crust change. In contrast, 
if the upper crust is relatively constant in thickness and the depth 
to the Moho is accomodated by varying the thickness of the lower 
crust, an inverse relationship is anticipated between gravity and mag- 
netic anomalies. 

Intra-layer lithologic variations caused by plutonic differentiates 
and regional and local metamorphism and alteration are known to be 
common in the upper crust and are almost certainly present on a wide- 
spread scale in the lower crust and upper mantle. Where ground-truth 
is available from direct observation of the upper crust, gravity and 
magnetic anomalies from these variations often are cor re la table . There 
are many exceptions, and even these are geologically informative, but 
in the main both the density and magnetization increase as the plutonic 
and metamorphic rocks become more mafic. Thus, we anticipate that 
lithogic variations within the lower crust will produce correlative 
anomalies. In fact, the correlation generally should be better from 
lower crustal rocks because slower cooling of the mafic rocks at depth 
will lead to increased exsolution of magnetite from the essentially 
non-magnetic components of the titanomagnetite series, tnus, increasing 
the magnetization of the denser mafic rocks. Furthermore, remanent 
magnetization which so conmonly complicates the correlation of upper 
crustal rocks, particularly the eruptive rocks, will be decreased in 
lower crustal rocks by the higher ambient temperatures and the increased 
age of the rocks. Thus, in general, intra-layer lithogic variations 
of the lower crust will yield directly correlative gravity and magnetic 
anomalies, although the correlations will not be universal or always 
direct. 

Finally, temperature variations, both regional and local within 
the crust, will affect the correlation of gravity and magnetic anomalies. 


Although the density of rocks will vary with temperature as a function 

of the volume coefficient of thermal expansion, the effect is minimal 

when considering the profound effect that temperature has upon magnetic 

properties. Magnetic susceptibility will be increased over certain 

temperature ranges and remanent magnetization will be decreased and 

viscous magnetization in the direction of the earth's magnetic field 

will be increased as temperature increases. However, the critical 

temperature is the Curie temperature, above which ferrimagnetic minerals 

are transformed into the essentially nonmagnetic, paramagnetic form. 

Thus, increased temperatures which reach the Curie point range of 
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ferrimagnetic minerals within the lower crust (500 -550 C) will 
decrease the thickness of this highly magnetic layer. The decreased 
magnetic anomaly will have an inverse relationship to the free-air 
gravity anomaly because regionally higher temperatures are associated 
with geodynamic processes which produce a thinned crust and an inflated 
elevation which may not be isostatically compensated. Both of these 
conditions will lead to an increase in mass and a relative positive 
free-air gravity anomaly correlating with a negative magnetic anomaly. 

Correlation of anomalies generally has been accomplished by sub- 
jective visual correlation of anomalies over regional areas. To 
decrease the subjectivity of the correlation techniques it is desirable 
to utilize a more quantitative approach, preferably a method based 
upon Poisson's theorem which theoretically relates magnetic and gravity 
anomalies from a conmon source. 

One such method is given by Chandler et_ aK , (1980) which employs 

Poisson '8 theorem in a moving-window linear regression analysis that 

traverses across multiple source gravity and magnetic data sets to 

quantitatively identify regions of spatially coincident anomalies. 

For regions of high correlation, the method also gives initial estimates 

of the physical property contrast ratios for the inferred sources. To 

investigate its applicability to the study of regional-scale potential field 

data, moving-window Poisson's analysis was performed between a number 

of gravity and magnetic anomaly maps prepared at satellite elevations 

o 

for North America and its marine environs south of 50 N latitude. 
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C. Quantitative Correlation Procedure 

Poisson's thsoram relates at an observation point the magnetic 
potential of a uniformly dense and magnetised body to the first 
derivative component of the gravity potential taken in the direction 
of magnetisation by the ratio of the body's magnetisation-to-density 
properties. Accordingly, Poisson's theorem can be extended to equate 
the total magnetic intensity anomaly reduced to normal (radial) polari- 
zation, AT(N), and the first radial derivative of the gravity anomaly, 
dAg/dr by the relation 

“»> ' * * 9 ? • " h,r * 

Ag » gravity anomaly, 

Am ■ density contrast of the source, 

C • universal gravitational constant, and 
Aj ■ magnetisation contrast of the source. 

Equation (IV-1) shows that the property contrast ratio (Aj/Am) for a 
single anomalous source can be obtained at a cotnnon observation point 
from its gravity and total magnetic field anomalies by calculating 
the radial gravity derivative and reducing the magnetic anomaly to 
radial polarization. 

For regional-scale , complex data sets containing interfering 
anomalies due to multiple sources, spatial segments may occur where 
the anomaly of a particular source is dominant. The approximate 
Poisson's relationship in such instances is given by the linear 
equation 

• A * etta * T? - (lv - 2> 

where the intercept A reflects the anomaly base level within the 
spatial segment considered, and the slope (Aj/Am) is a weighted 
average for the interfering sources that is most heavily weighted 
towards the dominant source. 

The concept of Poisson's analysis restricted to limited segments 
of a complex anomaly has been used successfully by Lyubimov and 
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Lyubimov (1968), Wilscn (1970), and Chandler et al., (1980). For com- 
bined analysis of multiple source gravity and magnetic data sets. 
Chandler et^al. , (1980) utilize the generalized Poisson's relationship 
(IV -2) in a moving-window linear regression analysis. The procedure 
is to center a window of specified size on each node of a uniform grid 
to which the appropriate gravity and magnetic data have been registered. 
Within the window, a least-squares linear regression is performed 
between the spatially coincident data where the radial gravity deri- 
vative serves arbitrarily as the independent regression variable. For 
each window position, the regression yields 1) a correlation coefficient 
that is indicative of the linear fit of the data within the window, 

2) a slope (Aj/Am) coefficient, and 3) an intercept coefficient (A) 
which gives the anomaly base level for the data within the window. 

The regression parameters are assigned to the central grid node of the 
window position. The window is shifted next one grid position and the 
regression analysis repeated. This procedure is continued until the 
entire data space has been subjected to moving-window Poisson's analysis. 

The size of the operational window chosen for a given appli- 
cation must be large enough that a minimum number of data points are 
involved to achieve a credible linear regression, yet small enough to 
resolve the anomalies of interest from the interference of neighboring 
anomalies. Currently, no quantitative method is evident for making 
this choice. However, Chandler et al. , (1980) have obtained useful 
results with windows that include at least 5 data points and are between 
0.5 and 1.0 times the ha If -width of the primary anomalies of interest. 

The behavior of the linear regression coefficients can be useful 
in establishing patterns of correlation between gravity and magnetic 
anomalies. The most significant linear regression exhibits correlation 
coefficients approaching either 1.0 or -1.0 which reflect, respectively, 
either direct (+) or inverse (-) correspondence between gravity and 
magnetic anomalies. Variations of the intercept parameter (A) reflect 
base level changes due to regional anomalies and the effects of 
localized anomaly interference. 
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For a particular window position* the slope of the regression 
line gives an estimate of (Aj/Am) for the dominant anomaly source. 
However* this ratio is not uniquely definitive of the anomaly source 
lithology because (Aj/Am) is a function of the physical properties 
of both the source and country rock. Also* different combinations 
of Aj and Am will give the same source ratio. The estimated slope 
value (Aj/Am) can serve* however* as an initial constraint for facil- 
itating further interpretation. 

Due to anomaly interference* variations of the slope values 
across multiple source data sets reflect both geologically signifi- 
cant and erroneous (Aj/Am) estimates. There is no quantitative 
criteria for discriminating geologically significant slope values. 
However, Chandler et al., (1980) recommend emphasis on slope values 
corresponding to the high-gradient limbs and peaks of spatially well 
matched gravity and magnetic anomalies which generally are character- 
ized by higher correlation coefficients and zones of spatially stable 
regression parameters. Particular emphasis should be placed on 
anomaly gradients because slight shifts of anomaly peaks will cause 
a marked decrease in the correlation over the peaks, but will not 
cause a marked change over the high-gradient areas. The dominance of 
a particular anomaly source is reflected in the spatial stability of 
the regression parameters. Spatial instability of these parameters 
generally characterizes transition zones between the dominance of 
adjacent anomalies. In such instances, the regression coefficients, 
particularly the intercept values, can exhibit marked local fluctua- 
tions. 

Spatial differentiation of equation (IV-2) shows that the moving- 
window Poisson's analysis readily can be extended to investigate 
correlations of the n-th magnetic spatial derivative anomaly field and 
the (n+1) -th gravity spatial derivative anomaly field. Spatial dif- 
ferentiation tends to emphasize the near-surface, higher frequency 
components of potential field anomalies. This can Le important to 
regional correlation analysis of satellite gravity and magnetic anomalies 
because spatial differentiation will minimize the interfering effects 
of gravity sources located deeper than the depth of the Curie isotherm 
which effectively defines the maximum depth extent for magnetic anomaly 


sources . 
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Suitable wavelength filtering in general can enhance the geologic 
significance of the results of moving-window Poisson's analysis by 
enhancing the spectral compatibility of gravity and magnetic anomaly 
data which contain the effects of common sources. Chandler et al., 
(1980) present a number of examples which illustrate the utility of 
wavelength filtering for moving-window Poisson's analysis. 


D. Data Set 8 and Processing 


To implement moving-window Poisson's analysis pertinent gravity 
anomaly derivatives and radially polarized magnetic anomalies and 
derivatives must be computed over a common observation grid. These 
processing procedures and the characteristics of the resultant data 
sets are considered below. 


1. Gravity Anomaly Maps 

The gravity maps used in this investigation all are derived from 
the smoothed, l°-averaged surface free-air gravity anomalies illus- 
trated in the contoured stereographic equal-area polar (SEAP) projection 
of Figure IV. 1. As described by Bowman et_ al. , (1979 ) , the smoothed 
gravity anomaly map was obtained by high-cut filtering anomalies with 
wavelengths less than about 8° from a data set of l°-averaged values 
provided by NASA-GSFC. The data set in Figure IV. 1 consists of (82,42) 

O O 

values uniformly spanning the region (219,5-300.5) E, (9.5-50.5) 14 
at 1° intervals. The amplitude range (AR) of the data set is between 
40 mgal and -62 mgal and the amplitude mean (AM) is -7.87 mgal. 

In general, the continent is characterized by relatively positive 
anomalies which trend roughly NS along the prevalent geological strike 
of the region. The Applachian and Ouachita Mountain Systems exhibit 
a slightly positive free-air gravity anomaly which continues south across 
Florida and Cuba. The Mississippi Embayment also is characterized by 
a weakly positive anomaly. The Cordilleran System is delineated by a well 
defined trend of positive free-air anomalies which attain maxima of 
roughly 30 mgal over Mexico and the Yellowstone geothermal region. 
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By contrast, oceanic anomalies are generally negative and trend 
EW along major fracture zones and trench systems. A number of fracture 
zones are well delineated by the EW striking components of the broad 
Pacific low (-20 mgal) off the west coast of California and the Baja 
Peninsula. The largest variations of the entire data set are the 
negative (-60 mgal) anomalies associated with the Puerto Rico Trench. 
The Gulf of Mexico also is characterized by a prominent negative 
(-50 mgal) free-air gravity anomaly. 

To project the smoothed free-air gravity anomalies to satellite 

elevations for comparison with POGO satellite magnetic anomalies, the 

data of Figure IV. 1 were related by least squares matrix inversion to 

a spherical grid of 435 gravity point poles. The grid of equivalent 

point sources was located at a depth of 400 km below the earth's 

o o 

surface and uniformly spanned the region (220-300) E, (10-50) N at a 
station interval of 2.857°. Computing the gravity effect of the 
point sources over the observation grid of Figure IV. 1 gives the 
equivalent source approximation of the smoothed free-air gravity 
"nomaly field shown in Figure IV. 2. A comparison of the amplitude 
means for the gravity anomaly maps of Figures IV. 1 and 2 shows that 
the equivalent point source approximation fits the original data with 
negligible error. 

To upward continue the smoothed free-air gravity anomalies to 
satellite elevations the equivalent point source field was computed 
next at 450 km elevation as shown in Figure IV. 3. Here the upward 
continued data set consirts of (81,41) values computed at integral 
degree nodes of the spherical grid uniformly spanning the region (220- 
300)°E, (10-50)°N at a station interval of 1°. Besides the obvious 
loss of resolution, Figure IV. 3 shows that the only positive anomalies 
surviving upward continuation are affiliated with the Cordilleran 
System of the U.S. , Mexico and Central America. Otherwise, the upward 
continued data are dominated by negative anomalies associated with 
the Puerto Rico Trench, the Gulf of Mexico, and the Pacific Plate. 

In Figure IV. 4, the first order radial derivative of the free- 
air gravity anomaly field at 450 km elevation is computed from the 
equivalent point source distribution for quantitative correlation via 
Poisson's theorem with satellite magnetic anomaly data reduced to the 
pole. The first radial derivative data provide enhanced resolution 
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of the higher frequency components of the upward continued gravity 
anomaly field. The Cordilleran System is characterized by a distinct 
positive zone of first order radial gradients which is flanked on the 
west by the broad negative zone of the Pacific Plate. To the east a 
well-defined zone of negative radial first derivative anomalies is 
evident which extends from the Great Lakes region SW across the mid- 
continent and terminates abruptly over the Gulf of Mexico as a large 
amplitude negative anomaly. Weakly positive first order gradients 
are affiliated with the Appalachian Mountain System which project 
SE across Florida and Cuba. The Atlantic is characterized by a broad 
zone of negative first order radial derivative anomalies with minimum 
values over the Puerto Rico Trench. 

The second radial derivative of the upward continued gravity 
anomalies is shown in Figure IV. 5 as derived from the equivalent point 
source distribution for comparison by a straightforward extension 
of Poisson's theorem with the first radial derivative of the satellite 
magnetic anomaly data reduced to the pole. The second derivative 
data provide further enhancement of the higher frequency (crustal ?) 
components of the upward continued gravity anomaly field. The 
spectral characteristics of the second derivative field are very 
similar to the characteristics of the surface free-air gravity anomalies 
of Figure IV. 1 which have been low-pass filtered for wavelengths of 
the order of 3° and larger. Accordingly, the second derivative 
anomalies show many of the associations previously noted for the 
smoothed l°-averaged surface free-air gravity anomalies of Figure IV. 1. 
In addition, the second radial derivative map of Figure IV. 5 shows 
a distinct arcurate banding of alternating positive and negative 
long-wavelength anomalies which extends westward across the Atlantic 
Ocean onto the eastern seaboard of the U.S. This arcuate banding of 
anomalies appears to project even further westward, although the 
pattern is disrupted frequently across the continental U.S. 

2. Magnetic Anomaly Maps 

POGO-1968 satellite magnetometer observations reduced to the pole 
at 450 km elevation are mapped in Figure IV. 6. This map was prepared 
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by least squares matrix inversion of multielevation satellite magnetic 
profiles on an equivalent sourcu distribution of spherical prismatic 
dipolar moments in the IGRF-1965 updated to 1968 according to the proce- 
dures described by Mayhew (1979) . Assuming these anomalies are derived 
primarily by induction, the satellite magnetic data were reduced dif- 
ferentially to radial polarization by recomputing the equivalent source 
fields using a uniformly polarizing amplitude of 60,000 gamma and normal 
geomagnetic field inclination at all source and observation points. 

The main features of Figure IV. 6 include the prominent EH trans- 
continental high which is breached in the region of the Mississippi 
Embayment. To the west, this high appears to overprint a NW trend of 
lows that corresponds to the Cordilleran System. The larger amplitudes 
of the transcontinental positive magnetic anomaly also correspond to a 
well-defined trend of gravity lows that extends from the Anadarko Basin 
of the eastern panhandle of Texas to the tri-state reyion of Indiana, 
Ohio and Kentucky (see Figures IV. 1 and 5) . This trend of gravity lows 
also characterizes a region of enhanced crustal thickness according to 
seismic evidence (Warren and Healy, 1973). 

Additional features of the radially polarized magnetic data of 
Figure IV. 6 include the high over the Anti lies -Bahamas Platform and 
the prominent magnetic low along the northern Gulf of Mexico which 
extends westward into Mexico. The magnetic relief of the Atlantic 
Ocean is considerable in comparison to the Pacific Ocean which appears 
as a region of relatively subdued magnetic character. 

Continental magnetic anomalies exhibit pronounced EH trends in 

contrast to the NS trends of the gravity anomaly fiolJ. This result 

is an apparent menifescation uf t»e longer wavelength components of 

the satellite magnetic data which are illustrated in Figure IV. 7. To 

construct Figure IV. 7, -he radially polarized data of Figure IV. 6 were 

o 

high-cut filtered tor wavelengths smaller than about 12 . The resultant 
map shows that major sources of positive long-wavelength magnetic 
anomalies include the previously discussed transcontinental region, the 
Antilles-Bahamas Platform and the southeastern coast of Mexico along 
the Middle America Trench. Negative long-wavelength magnetic anomalies 
characterize northern Mexico and the northwestern third of the U.S* 
where the western component of the latter anomaly is centered over the 
Yellowstone geothermal region. 
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For comparison with the second derivative map (Figure IV. 5) of the 
upward continued gravity anomalies# the first order radial derivative 
of the radially polarized magnetic anomaly field as shown in Figure IV- 
.8 is computed from the equivalent source field. To make the first 
derivative magnetic anomalies more compatible spectrally with the 
second derivative gravity anomalies# and thereby more amenable to 
correlation analysis via Poisson's theorem# the data of Figure IV. 8 
were high-cut filtered for anomaly wavelengths smaller than about 8°. 

The resultant smoothed first order radial derivative magnetic anomalies 
are illustrated in Figure IV. 9. 

The magnetic derivative maps provide an increased resolution of 
the previously mentioned anomalies. In addition, more subdued anomaly 
trends become apparent. For example# a series of negative anomalies 
is observed to extend along the continental side of the eastern sea- 
board of t.ne U.S. to Georgia where the trend turns westward along the 
southern margin of the U.S. and then across northern Mexico to the Gulf 
of California. On the oceanic side of the eastern U.S. seaboard, a 
trend of positive derivative anomalies is evident which projects from 
the Baltimore Canyon area southward until it intersects the positive 
anomalies of the Antilles-Bahamas Platform. Further eastward across 
the Atlantic Ocean, an arcuate pattern of banded negative and positive 
anomalies begins to emerge which shows a direct correspondence to the 
previously discussed gravity anomaly pattern of the Atlantic. The 
Pacific coast of North America by contrast generally is characterized 
by a trend of positive first order magnetic derivative anomalies. 

E . Quantitative Correlation Results and Discussion 

To investigate quantitative regional correlations of satellite- 
elevation gravity and magnetic anomalies over North America, moving- 
window Poisson's analysis was performed on the data sets described in 
the previous section. The results of these correlation analyses and 
their potential geological significance are outlined below. 

1. Correlation Analysis (1-5) 

Correlation analysis (1-5) was performed between the first radial 
gravity derivative map of Figure IV. 4 and the radially polarized magnetic 
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anomaly map of Figure IV. 6 using a 5x5 data point (»4°x4°) operational 
window. The positive and negative correlation coefficients derived 
from the moving-window Poisson's analysis are plotted separately in 
Figures IV. 10 and 11, respectively. These correlation coefficients 
jraphically characterize the oceans by direct correspondences of 
gravity and magnetic anomalies, whereas the continental correlations 
are predominantly inverse. 

A histogram of the correlation coefficients is given in Figure 
IV. 12 where the oceanic and continental anomaly correspondences are 
reflected as two distinct peaks of relatively high statistical cor- 
relation. Roughly 32% of the correlation coefficients have absolute 
values greater than or equal to ± 0.7, thus, indicating that a signi- 
ficant portion of the gravity and magnetic anomalies exhibit strong 
statistical correlation. 

An analysis of the correlation coefficient patterns in Figures 
IV. 10 and 11 shows that broad-scale regions of high correlation coef- 
ficients generally are associated with anomaly gradients. Correla- 
tions of anomaly peaks normally exhibit slight phase shifts of as much 
as a degree or two and, hence, tend to be characterized by fluctuating 
correlation coefficients. The slope and intercept coefficients as 
shown in Figures IV. 13 and 14, respectively, also have much the same 
pattern as the correlation coefficients. These results suggest that 
the regression parameters generally are more sensitive to phase shifts 
of the anomaly peaks than to phase shifts of the gradients. 

Accordingly, it is difficult to make broad geological generali- 
zations from these regression parameters except on an anomaly-by- 
anomaly basis. As an example, consider the previously discussed trans- 
continental region of negative gravity and positive magnetic anomalies. 
Much of this region is characterized by relatively stable intercept 
values and correlation coefficients which range between -0.6 to -0.8. 
The corresponding distribution of slope values, unich gives an initial 
constraint on the physical properties of the possible source for sub- 
sequent interpretation, is roughly equal to -0.06 oersted*cm 3 /gm. If 
the source of this gravity low and magnetic high corresponds to a 
thickened crust as suggested by previous arguments and seismic evidence 
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Figure IV. 12 Histogram of correlation coefficients derived from corre- 
lation analysis (1-5) . 
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(Warren « id Healy, 1973) , and the density contrast at the base of 

3 

the crust is -0.3 qm/ ’em , then this slope value indicates a magnetic 

-2 3 

polarisation for the source of roughly 1.8x10 emu/cm . This polari- 
sation value suggests source magnetisation that is several orders of 
magnitude greater than the value “ 5x10 emu/cm^ conmonly accepted 
for the lower crust (Rail, 1974) Shuey at al. , 1973). 

2. Correlation Analysis (II-5) 

Correlation analysis (IX-5) was performed between the second 
radial gravity derivative map of Figure IV. 5 and the first radial 
magnetic derivative anomalies of Figure IV. 8 using a 5x5 data 
point moving-window. This analysis should emphasise correspondences 
between nearer surface, crustal anomalies because both data sets are 
radial derivative maps. 

Positive and negative correlation coefficients derived from 
this analysis are plotted separately in Figures XV. 15 and 16. 
respectively. Again, these coefficients statistically characterise 
prevailing inverse and direct anomaly correspondences, respectively, 
for continental and oceanic regions. The histogram of the correlation 
coefficients given in Figure IV. 17 shows two broad peaks for the con- 
tinental and oceanic anomaly correlations. Approximately 23% of the 
anomaly correlations exhibit correlation coefficients with absolute 
magnitude greater than or equal to ± 0.7. However, inspection of 
Figures IV. 5 and 8 shows that a greater number of spatial correlations 
of anomalies appear evident than for the data sets considered in the 
previously described analysis. Hence, the decrease in the number of 
high correlations with respect to the results of analysis (1-5) sug- 
gests that the 4°x4° window may be too large to effectively resolve 
many of the higher frequency derivative anomaly correlations. 

Correlation coefficient fluctuations are similar to those 
observed for the slope and intercept values as shown in Figures XV. 18 
and 19. However, broad-scale variations of these coefficients 
generally correspond to anomaly gradients which make geologically 
significant generalizations unclear except on an individual anomaly 
basis. The results for the transcontinental anomaly correspondences 
are essentially the same as those obtained from the previous analysis. 
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Figure IV. 16 Negative (-1.0 to 0.0) correlation coefficients derived fro. correlation analysis IU-S) 
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figure IV. 18 Slope (Aj/^m) coefficients derived from correlation analysis (II-5) . 
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3. Correlation Analysis (III-3) 

Correlation analysis (III-3) using a 3x3 point (* 2°x2°) moving 
window was performed between the second radial gravity derivative map 
of Figure IV. 5 and the first radial magnetic derivative anomalies of 
Figure IV. 9 low-pass filtered for wavelengths roughly 8° and larger. 

The data sets considered here should be more spectrally compatible 
than those used in the previous analysis because the second radial 
gravity derivative map exhibits similar wavelength characteristics 
to the surface data of Figure IV. 1 which has been high-cut filtered 
for wavelengths roughly smaller than 8°. Also, the smaller operational 
window can be expected to better resolve the higher frequency cor- 
relations of the derivative anomalies. 

A number ot the derivative anomaly correlations are shown in 
Figure IV . 20 where selected second derivative gravity anomalies 
(heavy contours) are plotted on the smoothed magnetic anomaly 
derivative map. The spatial coincidence of the anomalies in Figure 
IV . 20 clearly demonstrates the continental area as a region chat- 
act erized basically by two types of inversely correlating gravity and 
magnetic anomalies. One type of inverse correspondence is exemplified 
by the positive gravity anomalies which are associated with breaches 
of the transcontinental magnetic high by the Cordillera and the northern 
Mississippi Embayment . The Yellowstone area also is characterized by 
well-defined positive gravity and negative maqnetic anomalies. The 
second type of continental anomaly correspondence is illustrated by 
the trend of negative gravity and positive magnetic anomalies which 
extends from the Anadarko Basin to the tri-state region of Indiana, 

Ohio and Kentucky. Oceanic anomaly correspondences are loss clear, 
although direct anomaly correlations generally appear to characterize 
the Gulf of Mexico, Caribbean and Atlantic Ocean regions of the eastern 
third of the study area. 

These generalizations, again, are supported statistically by the 
distributions of |>ositive and negative correlation coefficients derived 
from analysis (III-3) which are mapped in Figures IV. 21 and 22 , 
respectively. The histogram of the correlation coefficients illustrated 
in Figure IV. 2 3 shows the continental and oceanic anomalies exhibit 
relatively strong statistical correspondences. Roughly 41% of the 










101 


GRAVITY AND MAGNETIC ANOMALY 
CORRELATION ANALYSIS (m-3) 
CORRELATION COEFFICIENTS 
Z-4S0lm 


OCEANIC 



CONTINENTAL 


7ZZL 


ci bi uv 4 2Bb s*e 400 4 &s &itf set ei: 

DISTRIBUTION or CORRtLRTION COLT f I C I LNTS - 


Figure IV. 23 Histogram of correlation coefficients derived from corre- 
lation analysis (III-3). 
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derivative anomalies show correlation coefficients of absolute mag- 
nitude greater than or equal to ± 0.7. The increased number of relatively 
strong correlations reflects both the decreased site of the operational 
window and the enhanced spectral compatibility of the data sets. 

Fluctuations of the correlation coefficients tend to emphasize 
the slight phase shifts of the correlating anomaly peaks that are 
evident in Figure IV. 20. This behavior also is illustrated by the 
slope and intercept values which are plotted in Figures IV. 24 and 
25, respectively. Hence, despite the increased number of high statis- 
tical correlations, the relationship between the regression parameters 
and potential field anomalies remains essentially complex and must be 
considered on an individual anomaly basis. 

In the complicated behavior of the regression coefficients, however, 
continental inverse anomaly correlations tend to be characterized in 
a fairly uniform manner. Gravity lavs corresponding to the peaks of 
the transcontinental magnetic high all continue to exhibit correlation 
coefficients of -0.8 to -0.6 and slope values of the order of -0.06 
oersted*cm /gm. Gravity maxima corresponding to magnetic lows, on the 
other hand, as observed for the Yellowstone area and the breaches of 
the transcontinental magnetic high by the Cordillera and the northern 
Mississippi Embayment tend to all be characterized by correlation co- 
efficients up to -0.8 and slope values ranging between zero and -0.03 
oersted*cm 3 /gm . 

To investigate the possibility that the regression parameters con- 
tain significant long-wavelength information, the correlation, slope 
and intercept coefficients of analysis (III-3) were low-pass filtered 
for wavelengths roughly greater than 8° Positive and negative low-pass 
filtered correlation coefficients are given in Figures IV. 26 and 27, 
respectively, whereas the smoothed slope and intercept coefficients are 
illustrated, respectively, in Figures IV. 28 ar.d 29. The regression para- 
meters all exhibit similar fluctuating behavior which generally emphasizes 
the slight spatial disparity between correlating anomaly peaks. However, 
these variations are considerably simplified by the smoothing process and, 
hence, facilitate a more generalized interpretation of the correlating 
gravity and magnetic anomalies than the previously discussed results. 




Fiqure IV. 27 High-cut (X < 8°) filtered negative (-1.0 to 0.0) correlation coefficients derived from 
correlation analysis (III-3) . 
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For example, oceanic regions are broadly characterized by nail 
slope values in contrast to the continents. Zn consideration of the 
large magnetizations previously inferred for major portions of the 
continent, this result perhaps reflects decreased polarisation of the 
oceanic crust relative to continental crust because of the marked effect 
of temp e rature and time on the remanent magnetization properties of 
crustal extrusive rocks. The smoothed regression parameters also 
indicate more readily than in previous oxamples the uniform relation* 
ship between continental inverse anomaly correlations and the regres- 
sion coefficients. The transcontinental positive magnetic peaks, for 
example, which correspond inversely with gravity lows are uniformly 
characterized in the slope map of Figure XV. 28 by (Aj/Am) values 
ranging between -0.06 and -0.09 oersted»cm 3 /gm. 

The general characteristics of continental inverse gravity and 
magnetic anomaly correlations are further exemplified by the results 
for the profile along 37°N between 238°E and 284°E illustrated in 
Figure IV. 30. The top set of profiles (Figure IV. 30. A) shows the first 
radial derivative gravity data (dotted profile) and radially polarized 
magnetic data taken from Figures IV. 4 and 6, respectively, along the 
37°N parallel of latitude. The next set of profiles (Figure IV.30.B) 
corresponds to the second radial gravity derivative data (dotted profile) 
and smoothed first radial derivative magnetic anomalies along 37°N 
taken from Figures IV. 5 and 9, respectively. The third profile set 
(Figure IV. 30.0 gives both low-pass filtered (dotted profile) and 
unfiltered correlation coefficients corresponding to analysis (III-3) , 
and the final set of profiles (Figure IV.30.D) shows the slope values 
along 37°N extracted from Figures IV. 24 and 28 (dotted profile). 

These profiles clearly demonstrate the strong inverse correlation 
of regional gravity and magnetic anomalies that generally typifies 
North America. For the anomalies east of 270°E which *ire related to 
the gravity minima and magnetic highs discussed previously, the large 
(Aj/Am) values suggest crustal magnetization which is appreciably higher 
than indicated by conventional arguments, so that the magnetic high may 
be due to an intra-layer magnetic polarization variation. 
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Figure IV. 30 Gravity and magnetic anomaly correspondence at satellite 
elevations along 37°N latitude between 238°E and 284°E 
longitude. A) First radial free-air gravity derivatives 
(dotted profile) and radially polarized magnetic anomalies 
from Figures IV. 4 and 6 respectively. B) Second radial 
free-air air gravity derivatives (dotted profile) and 
smoothed first radially polarized magnetic anomaly deri- 
vatives from Figures IV. 5 and 9, respectively. C) Unfil- 
tered and high-cut filtered (dotted profile) correlation 
coefficients from correlation analysis (III-3). D) Unfil- 
tered and high-cut filtered (dotted profile) slope values 
from Figures IV. 24 anl 28, respectively. 
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The Mississippi Qnbayment is associated with a trail-defined 
gravity maximum and magnetic minimum. As described in section VI, 
seismic and surface gravity data indicate that the crustal structure 
of the Embayment consists of essentially four components. The major 
gravity source for the Embayment, however, corresponds to a portion 
of the lower crust with density contrast equal to 0.20 gm/cm . This 
particular feature also can account for the POGO satellite magnetic 

3 

data of the region if a polarization contrast equal to -0.0024 emu/cm 
is assumed. The resultant physical property contrast ratio (=-0.012 
oersted*cm 3 /gm) of the major gravity and magnetic source inferred for 
the Embayment compares very favorably with the slope value (<=-0.011 
oersted*cm 3 /gm) given in Figure IV.30.D. 

An especially broad band distribution of relatively strong inverse 

correlations is evident between the Mississippi Embayment and the 

Colorado Plateau. Uniform apparent (Aj/Am) values roughly equal to 

-0.037 oersted *cm 3 /gm tend to characterize this region, except for 

the marked local fluctuation in the vicinity of the 1° phase shift 

between the gravity anomaly peak at 259°E and the magnetic anomaly 
o 

peak at 260 E. 

The pronounced gravity maximum and magnetic low over the eastern 
margin of the Colorado Plateau probably reflects the increased heat 
flow and decreased crustal thickness of the region (Keller et al . , 1979} • 

Crustal thickening and lower heat flow over the Plateau result in a 
substantial decrease of gravity anomalies and slightly enhanced mag-* 
netics. However, west of the Plateau the effect of high heat flow 
and complicated tectonic styles yields an overall deterioration of 
anomaly correlations that is reflected in the general behavior of the 
correlation coefficients for this region. 

F* Conclusions 

Anomaly correlations as well as non-correlations provide useful 
information for the geologic interpretation of anomalous gravity and 
magnetic sources within the subsurface. To quantitatively assess 
patterns of anomaly correlation for multiple source data sets, Poisson's 

j 

1 

* i 
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theorem can be used to perform a regression analysis on appropriate 
gravity and magnetic anomaly data within a window that traverses the 
data sets. The regression parameters derived from this moving-window 
Poisson’s analysis include 1) correlation coefficients which provide 
a statistical measure of the coincidence and uniformity of the anomaly 
sources, 2) slope coefficients (Aj/Am) which provide initial con- 
straints on source physical properties to facilitate subsequent modeling 
for more detailed interpretation, and 3) intercept coefficients which 
reflect base level changes due to interfering anomalies. 

To investigate quantitative regional correlations of lonq-wave- 
length gravity and magnetic anomalies over North America and adjacent 
marine areas a number of data sets were prepared at 450 km elevation 
which satisfy the theoretical considerations of Poisson’s theorem in 
regions of common potential field sources. These data sets indicate 
the general relationship between gravity and magnetic anomalies over 
the continent is inverse. Regions of large relative crustal thickness 
such as the midcontinent appear to be characterized by negative gravity 
and positive magnetic anomalies at satellite elevations. Accordingly, 
this continental type of inverse anomaly correlation may reflect dis- 
placement of denser and less magnetic mantle material by an anomalously 
thick crust. However, the magnetic maxima probably also reflect an 
increase in the magnetic polarization of the crust in these areas. 
Negative magnetic and positive gravity anomaly correlations, on the 
other hand, as exemplified by the Yellowstone geothermal region may be 
related to regions of increased heat flow where the Curie point iso- 
therm is likely to be deflected to shallower depths to yield anomalous 
sources of net negative magnetization. 

The dis titict inverse characterization of continental gravity and 
magnetic anomaly correlations is supported statistically by the appli- 
cation of moving-window Poisson’s analysis to the regional data sets. 
These results indicate that, satellite elevation gravity and magnetic 
anomalies derived from continental crust demonstrate remarkable regional 
correlations that are viseful in decreasing the ambiguity of lithologic 
and structural interpretation. 
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Oceanic regions tend to exhibit overall poor anomaly correlations 
perhaps as a result of decreased crustal magnetization . Statistically, 
however, oceanic anomalies appear to be characterized by generally positive 
correlations which probably are best exemplified by the anomalies of 
the Caribbean, eastern Gulf of Mexico and the Atlantic Ocean. 

Poisson's moving-window regression analysis has been shown in 
this investigation to be particularly sensitive to phase shifts of 
anomaly peaks which complicates the interpretation of the quantitive 
correlations. Hence, in the application of this method it is desirable 
to obtain a high degree of spectral compatibility between the data sets 
to be correlated. This objective perhaps can be best achieved by 
simply filtering one of the data sets in the frequency domain using the 
spectral characteristics of the other data set as the filter. 

A number of reasons can be suggested which may account for the 
apparent phase shifting of anomaly peaks. Geologically significant 
mechanisms include sources with appreciable moments of remanent magneti- 
zation, sources with inhomogeneous properties of magnetization or 
density, and the differential response in time of magnetic and density 
properties to dynamic heat sources. Data reduction and processing 
problems also may lead to inadequate definition of anomalies. In the 
magnetic data set, however, this situation will soon be alleviated by 
the availability of scalar and vector magnetic anomaly data from the 
current Magsat program (Langel, 1979) which should provide higher 
resolution data to upgrade the quantitative correlation analysis of 
regional North American gravity and magnetic anomalies. 
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V. SPHERICAL EARTH GRAVITY AMD MAGNETIC ANOMALY MODELING 
BY GAUSS-LEGENDRE QUADRATURE INTEGRATION 




A. Abstract 


Increasing interest in the geologic significance of long-wave- 
length geophysical anomalies has focused attention on spherical gravity 
and magnetic modeling. Gauss-Legendre quadrature integration provides 
an especially powerful approach to the spherical earth modeling 
problem. A procedure has been developed using this numerical inte- 
gration scheme for calculating the anomalous potential of gravity and 
magnetic fields and their spatial derivatives on a spherical earth 
for an arbitrary body represented by an equivalent point source dis- 
tribution of gravity poles or magnetic dipoles. The distribution of 
equivalent point sources is determined directly from the coordinate 
limits of the source volume. Variable integration limits for an 
arbitrarily shaped body are derived from interpolation of points which 
approximate the body's surface envelope. The versatility of the 
method is enhanced by the ability to treat physical property varia- 
tions within the source volume and to consider variable magnetic fields 
over the source and observation surface. The method is particularly 
well suited for satellite gravity and magnetic anomaly modeling because 
the accuracy and efficiency of the application improves with increasing 
distance between source and observation point. A number of examples 
verify and illustrate the capabilities of the technique. 

B. Introduction 

Computation of theoretical anomalous gravity and magnetic fields 
due to geologic models is an important element in interpreting poten- 
tial field data and designing surveys. For regional gravity and mag- 
netic surveys measured in degrees of latitude and longitude, procedures 
are required which model directly in spherical coordinates potential 
field anomalies due to large-scale, arbitrarily shaped sources with 
variable density or magnetization characteristics. 
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In a review of computer modeling techniques, Bhattacharyya <1978) 
presents methods that in principle are suitable for spherical earth 
modeling of regional features. These procedures generally are based 
on approximations of the anomalous source as a group of prisms or 
polygonal laminae, the effects of which are evaluated and summed at 
each observation point to yield the total anomaly. However, for 
typical spherical earth modeling applications^ the bookkeeping problem 
involved with subdividing the large-scale source into simple forms to 
reflect arbitrary characteristics of geometry and physical properties 
commonly is formidable. 

A simple and more efficient procedure is to represent quadratures 
of the source volume by equivalent point sources according to the 
well known technique of Gaussian quadrature integration. In principle, 
the appropriate geometric distribution of equivalent point sources 
can be determined directly from the coordinate limits of the source 
volume, so that an accurate estimate of the source affiliated anomaly 
is obtained by evaluating and summing at each observation point the 
differential anomaly values due to each point source of the equivalent 
source distribution. 

Gaussian quadrature is a time-honored technique for numerical inte- 
gration that has been well studied in the general literature of numerical 
methods (e.g., Carnahan, 1969). Potential field modeling by Gauss- 
Legendre quadrature was used by Ku (1977) to evaluate gravity and 
magnetic anomalies in Cartesian coordinates due to bodies of arbitrary 
shape and magnetic polarization. In this discussion, the method is 
extended to spherical coordinates and the general problem of sources 
with arbitrary shape and variable density and magnetization properties. 

C. Description of the Method 

As illustrated in Figure V.l which outlines the main results of 
this section, it is convenient for gravity and magnetic modeling pro- 
blems to consider the anomalous body as being composed of a source 
volume distribution of gravity point poles or magnetic point dipoles, 
respectively. Hence, to estimate the anomalous gravity or magnetic 
field at some observation point, it is necessary to evaluate and sum 
the anomaly values due to each of these point sources at the observation 
point. 
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Figure V.l Gravity and magnetic anomaly modeling in spherical coor 
dinates cf a geologic body with arbitrary shape and 
physical properties by Gauss-Legendre quadrature inte- 
gration (see text for details) . 
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In particular, the radial anomaly, Ag, due to a gravity point 
pole referenced to a geocentric coordinate system as shown in Figure 
V.l is given by 

d 1 an 

Ag ■ (-)^p)*Am, where (V-l) 

.o 3 2 

6 ■ universal gravitational constant (>6.67*10 cm /(gm*sec )), 

R ■ distance from observation point (r, 0 ,$) to source point (r' , 6 ' ,$ ' ) , 

r ,r ' - radial distances from earth's center to the observation point 
and source point, respectively, 

6 ,8' ■ co-latitudo coordinates of observation and source points, 
respectively, 

♦ >♦' > longitude coordinates of observation and source points, 
respectively, and 

Am = mass contrast of the point pole. 

For a magnetic point dipole, on the other hand, the total magnetic 
intensity anomaly, AT, is given by 

AT > {-u-t(u'-$'(i))}*Aj, where (V-2) 

* gradient operators in observation point and source point 

coordinates, respectively, j 

i 

Aj » magnetization contrast of the point dipole, and j 

•f + 

u,u' = unit geomagnetic field vectors at the observation point and j 

source point, respectively. I 

Conventionally, the unit vectors u and u' are expressed in terms of 

geomagnetic field inclination (1,1') and declination (D,D'). Also, 

. i 

when the point dipole polarization is by induction 

Aj ■ Ak*F', where (V-3) 

Ak = magnetic susceptibility contrast of the point dipole, and 
F' > scalar geomagnetic field intensity at the source point. 

Hence, in regional-scale applications of equation (V-2) geomagnetic • j 

l 

field models such as the XGRF-1965 (Cain et al . , 1967) normally are $ 

used to obtain pertinent values of (I',D',F') at the source point 
and (I ,D) at the observation point. 

Consideration of equations (V-l) and (V-2) shows that the incre- 
mental gravity and magnetic anomaly values due to a point source located 
a distance R from the observation point can be generalized according 
to the relation 
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q(R) *Ax, where (V-4) 

q(R) ■ the geometrical point source function which essentially 
describes the inverse distance between the observation 
point and source point, and 

Ax ■ the appropriate physical property contrast of the point source. 
Hence, to determine the total potential field anomaly it is necessary 
to evaluate at each observation point (r ,6 ,♦) the volume integral given 
by 

Xtb $> (q(r,0,$ jr' ,0' ,$')*Ax)dr'd0 *d$ ' , (V-5) 

*a a a 

where tne primed variables refer to the coordinates of the anomalous 
body such that 

' - lower and upper longitude limits of the source volume, 

A D 

0’,0* • lower and upper co-latitude limits of the source volume, end 
a o 

r',r' * lower and upper radial limits of the source volume, 
a o 

Consider, now, the numerical evaluation of the innermost integral 
in equation (V-5) . Most numerical integration techniques involve the 
use of interpolation polynomials to approximate the integrand according 
to a summation formula of the general type given by 

r . ni 

£,b (q(r')*Ax)dr' * £ A^qtrp^Ax, (V-6) 

a i«l 

where the ni values A^ are the weights to be given to the ni func- 
tional values (q(r£)*Ax) evaluated at the interpolation coordinates, 
r^. Conventional integration formulae such as Simpson' s rule or the 
trapezoidal rule assume equal spacing of the arguments r£ which is 
generally appropriate when dealing with an integrand that is not well 
known analytically. However, as shown in equations (V-l) and (V-2) 
the integrand being considered here involves a familar analytic function 
which may be computed for any argument to great precision. In such 
instances, Gaussian quadrature formulae can be developed to yield 
selected values of interpolation rj^ and coefficients A^ so that the 
sum in equation (V-6) gives the integral exactly when (q(r')*Ax) is a 
polynomial of degree 2ni or less (e.g., Carnahan, 1969). 

In general, it may be shown that the Gaussian coefficients A^ can 
be obtained from a polynomial of order ni which is orthogonal over 
the interval of integration such that the ni points of interpolation 
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r| are the zeros of the polynomial. Families of orthogonal polynomials 

which commonly are used to develop Gaussian quadrature formulae include 

Legendre# Laguerre, Chebyshev and Hermite polynomials. However# in 

this discussion only Legendre polynomials are considered which is 

the prototype of the Gaussian method. 

Legendre polynomials P , (r') of order ni which are orthogonal 
** ni 

over the interval -1 < r' <1 are given by 

P . (r') - (-—7 ) (— (r ,2 -l) ni ), Where P <r») - 1. (V-7) 

ni 2 ni niJ d;* ni ° 


Accordingly, the standard Gauss-Legendre quadrature over the interval 
(•1,1) is given by 




ni 

(q(r')*Ax)dr' * E A *q(r')*Ax, 
i-1 1 


(V-8) 


where the interpolation points r| at which the integrand is evaluated 
are the zeros of equation (V-7) and the Gaussian coefficients are 


A 


i 


2 d-?' 2 ) 

^ P ni-^i )}2 ' 


(V-9) 


Now for arbitrary limits of integration such as indicated in 
equation (V-6) it is necessary to map the standard interval -1 <r^ <1 
into the interval of integration r' <rj <r' according to the trans- 
formation 


r! (r'-r' ) + (r'+r‘) 
r . » k b ■ *> a. 

i 2 


(V-10) 


Accordingly , the integral in equation (V-6) can be approximated as 


£fb (q(r‘ ) *Ax)dr ' ■ — (q(~ — - — - 


r'(r'-r') + (r>r‘) 


b a 


■}*Ax)dr' 


(r'-r* ) ni 

« — — — E A.*q(r!)*Ax. 
i-1 1 


(V-ll) 


Extending this procedure to volume integrals is straightforward. Thus# 
the Gauss-Legendre formula for the general evaluation of equation (V-5) 
is given by 
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fib fib fib fq(r' ,8* # +')*Ax)dr 9 d6'd$ 9 
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where 


n 

9 i 
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r£,8 9 ,$ 9 - coordinates of the subdivision in the limits of integration 
3 from -1 to 1 which correspond to zero nodes of equation CV-7), 

AftAjtAk ■ Gauss-Legendre quadrature coefficients given by equation (V— 9 ) , 

$ ' ' ■ lower and upper volume limits of the body in the k-th 

a coordinate of longitude* 

6 9 *8 9 * lower and upper volume limits of the body in the j-th 

^ a ^ coordinate of co-latitude* and 

r ia ,r ib " lower aru * upper volume limits of the body in the i-th 
a radial coordinate. 

The quadrature formula given in equation (V-12) shows that gravity 
and magnetic anomalies can be accurately computed by summing at 
each observation point the anomalous effect of nk*nj*ni equivalent 
point sources located at source point coordinates (r| ,0 ^ . where 
each of the differential point source anomalies is appropriately 
weighted by Gauss-Legendre quadrature coefficients and the volume coor- 
dinate limits of the anomalous body being modeled. This result is 
well suited ior machine computation where the input consists prin- 
cipally of values of the integrand for selected source points (r^,0^ , 
affiliated subdivision coordinates (r^*0 * and coefficients 
(A^*A^«A ]( )* and the volume coordinate limits of the body for each 
dimension of every source point coordinate (r|, 8 *,$£). 

The selected values of the integrand for gravity or magnetic 
modeling purposes are readily obtained from equations (V-l) or (V-2)* 
respectively. Also, the Legendre subdivision coordinates of the interval 
(-1,1) and associated Gaussian coefficients may be computed directly 
from equations (V-7) and (V-9) , respectively; although* it generally is 
found to be more machine efficient to input these values from tables 
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using an algorithm such as dascribad by Carnahan (1969) . Values which 
are applicable for such algorithms are tabulated to 30 digit precision 
for orders n ■ 2 to 512 by Stroud and Secreat (1966). However , for 
most regional lithospheric modeling applications experience suggests 
subdivision coordinates and associated coefficients to 10 digit pre- 
cision for orders up to n ■ 16 normally are sufficient. 

For a uniformly dimensioned body such as a prism, the integration 

limits for the evaluation of equation (V-12) are easy to specify. Zn 

this case, for example, “ (9^#9£) # * *®a ,0 b* ftnd 

(r! ,r' ) » (r',r'). However, for the more general case of a body with 
la lfi a o 

arbitrary shape the integration limits are known in one dimension only 
so that the problem of determining the integration limits in the remaining 
two dimensions for each equivalent point source coordinate must be 
considered. 

Procedures can be developed to efficiently handle this problem, 
such as the method described by Ku (1977) where a modified cubic spline 
function is used to interpolate the desired integration limits from a 
set of body point coordinates which provide a rough approximation of 
the surface envelope of the body. Typically, the procedure is to 
specify, for example , the longitudinal limits of integration of the 
body to obtain the nk Gauss-Legendre nodes as described above. 
Interpolations of the body point coordinates are performed next to 
determine the maximum and minimum latitude coordinates of the body 
for each longitude coordinate $£. These values of course provide the 
latitude limits of integration for evaluating the nj nodes 6^. Simi- 
larily, the radial coordinates of the body points are interpolated at 
each horizontal coordinate ($£,0 to yield appropriate radial limits 
of integration from which the ni nodes can be determined. Pro- 
cedures such as this are readily adapted for efficient machine processing 
so that the integration limits of arbitrarily shaped bodies can 
be determined accurately for evaluation of the quadrature formula 
given in equation (V-12) . 

The quadrature formula in general has considerable versatility 
in modeling applications because anomalous gravity and magnetic 
potentials and their respective spatial derivatives of any order are 
all linearly related. Hence, to model the radial derivative of the 
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potential field anomaly due to an arbitrary source, for example, it 
is necessary simply to exchange q(R) for 3q(R)/dr in equation (V-12). 
Additional geophysically interesting quantities which can be modeled 
from simple linear transformations of the Integrand of equation (V-12 ) 
include the anomalous potential, vector anomaly components and the 
spatial derivatives of any order. Relative geoidal anomalies can be 
modeled as well by computing the anomalous gravitational potential 
of the body at the surface of the earth and dividing it by normal 
gravity according to Brun's formula (Heiskanen and Moritz, 1967). 

Bodies with variable physical property contrast also are accom- 
modated readily by the quadrature solution. To emphasize the fact that 
Ax can be expressed as a function of source position (i.e., Ax ■ 

Ax(r' ,9' ,♦')) when necessary, the generalized point source anomaly 
(q(R) *Ax) has been carried intact throughout the foregoing developments. 
Conventional modeling procedures normally assume bodies with uniform 
physical properties so that each property variation must be modeled 
as a separate body. To model the total magnetic intensity anomaly due 
to a regional source subject to regional variations of geomagnetic 
field induction by the well known method of Talwani (196S) , for example, 
requires that the source be subdivided into blocks wherein the geo- 
magnetic field polarization is uniform. Bach block in turn must be 
subdivided into a group of approximating polygonal laminae to insure 
accuracy of the numerical integration which evaluates the anomaly. By 
contrast , the Gauss-Legendre quadrature approach is to simply polarize 
each point dipole in the quadrature formula according to the polari- 
zation characteristics of an acceptable numerical model of the gexnag- 
netic field (e.g., the IGRF-1965 by Cain et al., 1967)) This e r pie 
also illustrates the technique for modeling bodies with remanent 
magnetization by quadrature formulation. Here, the magnetization of 
the point dipoles is achieved using a polarization field that represents 
the vector sum of induced and remanent magnetic polarizations. 

Finally, considerable computational flexibility is available for 
the practical implementation of the quadrature formulation Ku (1977) 
noted the tradeoffs which occur in applications of the raetr between 
efficient computation and the accuracy of the solution. A large 
number of Gauss-Legendre nodes insures a very accurate quadrature 
solution, although the same degree of accuracy often can be achieved 
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by a substantially smaller number of equivalent point sources. In 
fact, the accuracy of the solution remains essentially unchanged for 
different numbers of nodes as long as the node spacing is smaller 
than the distance to the observation point. Hence, the minimum number 
ct aodes to specify in a given application should be such that the 
distance to the observation point is greater than the node spacing. 

In general, then, the accuracy of the quadrature formulation can readily 
be controlled by adjusting the number of Gauss-Legendre nodes relative 
to the elevation of the observation point. 

D . Illustration the Method 

A computer program is described in von Frese et_ al. , (1980) which 
was developed for regional lithospheric gravity and magnetic modeling 
applications by Gauss-Legendre quadrature. The program was used to 
construct the examples of this section which illustrate some of the 
capabilities of the method. 

To demonstrate and verify the method for regional scale modeling 
applications, gravity and magnetic anomalies due to a three-dimensional 
spherical prism were modeled n a spatial scale small enough that the 
results could be compared readily with conventional modeling techniques 
in Cartesian coordinates. In Figure V.2, for example, a comparison of 
gravity anomalies is given due to a 6.67 km thick prism with density 
contrast 0.25 gm/cm^. The observation grid consists of (16,16) stations 
uniformly spanning the region (273-274) °E, (40-41) °N at an elevation 
(Z) of 6.67 km above the top of the prism. 

The well known method of Talwani and Ewing (1960) was used to 
calculate the gravity effect of the prismatic model as shown in Figure 
V.2. A. The gravity anomaly was determined by evaluating and summing at 
each observation point the gravity anomalies due to 11 horizontal poly- 
gonal laminae used to approximate the prism. As shown in Figure V.2. A, 
the resultant anomaly has an amplitude range (AR) between 36.0 mgal 

and 0.3 mgal and an amplitude mean (AM) of 5.4 mgal. The anomaly 

o 

was computed and contoured in Cartesian coordinates assuming 1 = 100 km. 

This assumption distorts slightly the true geometry of a spherical 
prism in this 1° x 1° region in northwestern Indiana where 1° in longi- 
tude and latitude is more nearly equal to 88 km and 112 km, respectively. 
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In Figure V.2.B, the gravity anomaly of the spherical prism confuted 
by Gauss-Legendre quadrature integration is contoured on a stereo- 
graphic equal-area polar (SEAP) projection. In this case, the gravity 
anomaly was calculated in spherical coordinates by evaluating at each 
observation point an nk*nj*ni « 8*8*2 « 128 point quadrature formula 
where the integration limits were specified directly from the spherical 
coordinate limits of the prism volume. The results shown in Figure V.2 
indicate the two methods agree very well with respect to the relatively 
small spatial scale under consideration. The slight disparity between 
the results probably is due to geometric distortion introduced by the use 
of Cartesian coordinates for calculation of the spherical prism anomaly. 

A comparison of magnetic anomaly fields is illustrated in Figure V.3 
for the prism model wi H h volume magnetic susceptibility contrast* equal 
to 0.0005 emu/cm' 3 and uniform polarizing intensity F* = 60,000 gamma, 
inclination I* = 75° and declination D* = 0°. Over the observation 
grid, uniform geomagnetic field attitude characteristics of inclination, 

I * 75° and declination D = 0° also are assumed. The total magnetic 
intensity anomaly of the spherical prism calculated in Cartesian coor- 
dinates according to the method of Talwani (1965) is shown in Figure V.3. A. 
The resultant anomaly corresponds well with the magnetic anomaly cal- 
culated in spherical coordinates by Gauss-Legendre quadrature in Figure 
V.3.B. This correspondence also holds for the exact closed form solution 
for the magnetic anomaly of the prism in Cartesian coordinates given by 
Andreasen and Zietz (1969) which is illustrated in Figure V.4.A. 

Figure V.4.B illustrates the magnetic anomaly calculated by Gauss- 
Legendre quadrature integration using the algorithm described previously 
to interpolate the limits of integration from a set of points which 
sample the surface envelope coordinates of an arbitrarily shaped body. 

Here, the set of body points used in the algorithm comprise the longitude, 
latitude, and top and bottom radial coordinates of the prism for each 
node of the observation grid that overlies the boundary or interior of 
the prism. The actual interpolation was performed using the "mean slope" 
cubic spline function described by Ku (1977) . The excellent agreement 
between Figures V.4.B and V.3.B shows that the algorithm can be used to 
obtain interpolated integration limits for arbitrarily shaped bodies with 
sufficient precision to insure accurate anomaly determinations. 


*) 1 emu/ cm 3 * l*4ir 





anomaly comparisons for a spherical prism where the Gauss— Legend re quadrature 
(B) is developed from interpolated limits of integration. 





3 < 


Some of the versatility of potential field modeling by the Gauss- 
Legendre quadrature formulation is illustrated in the next two figures. 

In Figure V. 5 , the radial and horizontal components of the magnetic anomaly 
are shown due to the spherical prism model. These anomaly components 
were obtained by evaluating and summing at each observation point the 
vertical and horizontal components of the point dipole anomalies 
(i.e., equation (V.2)) due to the 128 equivalent point sources used in 
the quadrature formula. Other geophysically interesting characteristics 
of the prism's anomalous gravity and magnetic potentials can be deter- 
mined in this fashion as well, including the associated geoidal anomaly, 
vector anomaly components and spatial derivatives of any order. 

In Figure V.6, the effect of including a remanent component in the 
polarization of the prismatic model is demonstrated. The polarization 
of a source represents the vector sum of its remanent magnetization 
field and the inducing geomagnetic field. For the case considered 
here, the effect of remanence yields a uniform polarization field 
with inclination I* = 30° and declination D* = N30°E in contrast to 
the attitude of the geomagnetic field over the observation grid which 
is I = 90° and D = 0°. The anomaly for this case in Cartesian coor- 
dinates according to the closed form solution given by Andreasen and 
Zietz (1969) is illustrated in Figure V.6. A. Comparing this result 
with the spherical coordinate solution given by Gauss-Legendre quadrature 
in Figure V.6.B again demonstrates good correspondence with respect to 
the small spatial scale being considered. 

The foregoing examples show that the formulation of the quad- 
rature solution in spherical coordinates for bodies of arbitrary 
shape and physical properties gives results which reduce to well studied 
solutions under relevant modeling circumstances. Accordingly, these 
examples help provide useful practical incentives for implementing 
potential anomaly field modeling of regional geologic sources by Gai^ss- 
Legendre quadrature integration. Ku (1977) also presents a number of 
comparisons, although in Cartesian coordinates, which further illus- 
trate the excellent agreement between the quadrature formulation and 
other published solutions for such idealized sources as the step model, 
rectangular prism, sphere, and vertical and inclined dykes. 
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As a mors regional-scale example, consider the application of Gauss- 
Legendre quadrature integration for modeling gravity and magnetic anomaly 
signatures at the satellite elevation of 450 km due to the crustal 
thickness anomaly which is gridded in Figure V.7.A. This zone of enhanced 
crustal thickness, extending from the Texas panhandle northeastward 
into Kansas, is portrayed by seismic evidence (Warren and Healy, 1973) 
as roughly 10 km of lighter crustal material displacing denser (mantle) 
material to a depth of 50 km. To illustrate the gravity modeling pro- 
cedure, a density contrast of -0.3 gm/cm^ was assumed for this feature. 

Magnetic considerations, on the other hand, suggest that the crustal 

thickness anomaly also may represent a zone of positive magnetization 

contrast due to downward deflection of the Curie isotherm in the region 

of thick, cooler crustal material or thicker magnetic :rust extending 

into non-magnetic mantle material. For the purposes of this exangple, 

-3 3 

a volume susceptibility contrast of 10x10 emu/cm was assumed for the 
zone of enhanced crustal thickness which is representative of the general 
magnetization reported for the lower crust (e.g.. Hall, IS 74? Shuey 
et al. , 1973) . 

To determine the variable limits of integration for the evaluation 
of the quadrature formula, the crustal thicKness anomaly was referenced 
to the body point grid shown in Figure V.7.A. For each Gauss-Legendre 
node in the <J>-longitude coordinates of the body, the 6-latitude limits 
of the body were interpolated from the 6 points describing the body's 
boundary. The radial limits, in turn, were interpolated from the 6 
boundary points and 3 interior points of the body for each Gauss- 
Legendre node with horizontal ($, 0) -coordinates. For the particular 
model considered here, of course, the radial coordinates of the 9 body 
points used to approximate the subsurface configuration of the body 
were specified to represent a uniform thickness of 10 km. 

The resultant gravity anomaly at 450 km elevation due to this 
feature is illustrated in Figure V.7.B. The gravity effect was cal- 
culated by evaluating an nk*nj*ni * 16*16*2 « 512 point Gauss-Legendre 
quadrature formula over the (41,41) grid of observations. The magnetic 
effect of the crustal thickness anomaly also was computed in this manner 
at 450 km elevation. The resultant total magnetic intensity anomaly in 
the IGRF-1965 updated to 1968 is demonstrated in Figure V.8.A. To 
remove anomaly distortion due to regional variability of the 
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(10 km) is indicated by the 50 km depth contour on which the grid is superimposed. B) 
SEAP projection of the gravity anomaly for the crustal thickness model calculated at 450 
elevation by Gauss-Legendre quadrature integration. 




CRUSTAL THICKNESS 



Relative precision of the gravity anomaly (A) and total magnetic intensity anomaly (B) 
calculations is essentially unchanged at 450 km elevation by a 93.75% reduction in the 
number of Gauss-Legendre nodes. 
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reference field, the magnetic effect wae computed next assuming a 
uniform polarising field intensity of 60,000 gamma and radial geo* 
magnetic field inclination at both source and observation points. 
Accordingly, Figure V.8.B illustrates the resultant magnetic anomaly 
field reduced to radial polarisation at 450 km elevation for the 
crustal thickness model. 

Finally, Ku (1977) noted that for accurate implementation of the 
method the distance between the equivalent source points and the obser- 
vation point must be greater than the equivalent source point spacing 
within the bodv. This limitation can be minimised in practice by either 
subdividing the body into an appropriate number of smaller bodies, 
or increasing the number of equivalent point sources, or increasing 
the distance between the observation point and the body. The latter 
consideration suggests that the Gauss-Legendre quadrature formulation 
is especially well suited for modeling satellite leval gravity and 
magnetic anomalies because of the large elevations involved with these 
measurements. 

The principle is demonstrated in Figure V.9 where the gravity 
anomaly and total magnetic intensity anomaly due to the crustal thick- 
ness model are computed at 450 km elevation using an nk*nj*ni » 4*4*2 ■ 
32 point quadrature formula rather than the 512 point formula. The 
comparison of the mean gravity anomaly amplitudes of Figures V.7.B and 
9. A shows accordingly that a 93.75% reduction in the number of equi- 
valent point sources results in only a 1.67% reduction in the relative 
precision of the gravity anomaly calculation. A slightly smaller 
reduction of 1.37% is evident for the relative precision of the magnetic 
anomaly calculations illustrated in Figures V.8.A and 9.B. 

B. Conclusions 

Regional gravity and magnetic anomaly modeling for arbitrarily 
shaped lithospheric sources with variable physical properties can be 
achieved accurately and efficiently in spherical coordinates using 
Gauss-Legendre quadrature integration. The procedure involves repre- 
sentation of the anomalous source as a distribution of equivalent point 
gravity poles or point magnetic dipoles which contribute incremental 
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anomaly values that are evaluated and stunned at each observation point 
to obtain the total anomaly. The distribution of equivalent point 
sources is determined directly from the volume limits of the anomalous 
body. For an arbitrarily shaped body, the variable limits of inte- 
gration can be obtained from interpolations performed on a set of 
body points which approximate the surface envelope of the anomalous 
source. 

A chief practical advantage of the method is its considerable 
versatility. The physical properties of the emiivalent point sources, 
for example, can be individually varied to reflect physical property 
variations of the body that is being modeled. The method also can be 
extended readily to model the geoidal anomaly, vector components, and 
spatial derivatives to any order of the body's anomalous gravitational 
and magnetic potentials. Finally, the accuracy of the method can be 
controlled by adjusting the number of equivalent point sources or 
the distance between the source and observation point. In this regard, 
the method is particularly well suited for satellite gravity and 
magnetic anomaly modeling because the efficiency and accuracy of the 
application increases with increasing distance between source and 
observation points. 

In consideration of the foregoing results, it is concluded that 
Gauss -Legendre quadrature integration facilitates a powerful and 
efficient approach to spherical earth modeling of regional-scale litho- 
spheric gravity and magnetic anomaly sources. Accordingly, the method 
has widespread application in the analysis and design of regional -scale 
gravity and magnetic surveys for lithospheric investigation. 
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VI. GRAVITY AND MAGNETIC ANOMALY MODELING OF MISSISSIPPI 
EMBAYMENT CRUSTAL STRUCTURE AT SATELLITE ELEVATIONS 
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A. jjjstract 

A model for the three-dimensional crustal structure of the northern 
Mississippi Embayment is generalized from published surface wave disper- 
sion, seismic refraction, and gravity studies of the region. The gravity 
and magnetic anomaly signatures due to this model are computed at 450 km 
elevation by Gauss-Legendre quadrature integration for comparison with 
observed anomalies at satellite elevations. The computed positive gravity 
anomaly compares well with upward continued free-air gravity data sug- 
gesting that the generalized model is representative of the crustal 
structure of the Embayment. Magnetic anomaly calculations show that the 
pronounced minimum observed over the Embayment in the POGO satellite 
magnetometer data can be accounted for by a decrease in the magnetization 
of the lower crust which corresponds to the major gravity source of the 
region. The results of this investigation support the failed-rift ^ 

hypothesis for the origin of the Mississippi Embayment. Accordingly, 
these results suggest that observable gravity and magnetic anomalies 
characterize failed rifts (aulacogens) at satellite elevations, where 
the primary source of both anomalies is a high density rift component 
of non-magnet ic lower crustal material. 

B. Introduction 

Satellite-level magnetic and gravity anomaly interpretational methods 
largely have been limited to visual spatial correlation of anomalies 
with known geological features and simple flat earth modeling approxi- 
mations. However, the full importance of satellite-level anomalies, 
which map the anomalies of large scale geologic features that are difficult 
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to recognize and isolate in low- level data, will not be achieved until 
direct modeling of large-scale structural features is accomplished 
utilizing the spherical earth condition. Modeling of geologic features 
incorporating available geologic and geophysical data as constraints 
can be employed as an interpretive technique by altering tie model until 
the computed anomalies match the observed anomalies. Furthermore, 
spectral analysis of the calculated satellite-level magnetic and gravity 
anomalies of geologic models provides critical information for designing 
filters to isolate the observed anomalies of these features. 

In the previous discussion a technique was described and demon- 
strated for numerically modeling potential field characteristics due to 
regional-scale geologic sources by Gauss-Legendre quadrature integration. 
Herein, the method is used to obtain a preliminary view of the gravity 
and magnetic anomaly characteristics for the Mississippi Embayment at 
450 km elevation. This information, in turn, is pertinent to evaluating 
the feasibility of using satellite gravity and magnetic surveys for 
detecting anomaly signatures due to failed rifts. The Mississippi Embay- 
ment is particularly suited for this application because a number of 
geophysical constraints are available for developing a first-order model 
of its crustal structure. 

C. Model Description 

The Mississippi Embayment represents a broad, spoon-shaped reentrant 
of Mesozoic and Cenozoic sedimentary rocks which extends into the 
Paleozoic terrane of the North American craton from the south. As 
shown in Figure VI. 1, the axis of this feature roughly parallels the 
Mississippi River tapering northward into the tectonically active region 
of the New Madrid seismic zone. 

Based on the regional geology, Burke and Dewey (1973) suggested 
the embayment is a Mesozoic aulacogen developed from a triple junction 
located in the vicinity of Jackson, Mississippi. However, an integrated 
analysis of gravity, seismic, stratigraphic and petrologic data by 
Ervin and McGinnis (1975) suggests the embayment is a late Precambrian 
aulacogen which was reactivated most recently in the late Cretaceous by 
tensional forces initiated during the formation of the present Atlantic 
Ocean basin by subsidence of the Gulf Coastal Plain. 
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Index map of the Mississippi Embayment (shaded contour) 
region where G is the gravity profile studied by Ervin 
and McGinnis (1975) , R is the seismic refraction line 
studied by McCamy and Meyer (1966) , and S is the surface 
wave propagation path studied by Austin and Keller (1979) 
The lateral subsurface boundary of the model generalized 
in Figure VI. 2. C for this study also is indicated. 
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Austin and Keller (1979). C) The density contrast model generalized for this study from 
Austin and Keller (1979). Body //2 with magnetization contrast -0.0024 emu/cm 3 also was 
used to model POGO satellite magnetometer data over the Embayment . D) The gravity data 
of Cordell (1977) used to develop the lateral subsurface configuration outlined in Figure 
VI. 1 of the density contrast model generalized in C. 




Figure VI. 2. A is a cross-section of the density structure of the 
Mississippi Enuuayment given by Ervin and McGinnis (1975) along a profile 
between Yellvillo, Arkansas and Scottsboro, Alabama (hereafter called 
Y-S profile). This density model was synthesised from regional gravity 
data derived from the U.S. Bouguer gravity anomaly map of Woollard and 
Joesting (1964) and the results of a reversed seismic refraction profile 
between Little Rock, Arkansas and Cape Girardeau, Missouri as described 
by McCamy and Meyer (1966) . Austin and Keller (1979) integrated the work 
of McCamy and Meyer (1966) with an analysis of Rayleigh wave dispersion 
along a propagation path between Oxford, Mississippi and Florissant, 
Missouri to obtain a similar density model for the Y-S profile which is 
illustrated in Figure VI. 2. B. An index map for locating these various 
studies is given in Figure VI. 1. In general, the crustal cross-sections 
shown in Figures VI. 2. A and B support the failed-rift model for the origin 
of the Mississippi Embayment. 

The agreement of surface wave, seismic refraction, and gravity data 
in the region of the Embayment suggests that the crustal cross-section 
given in Figure VI. 2. B can be useful for developing a reasonably valid 
three-dimensional model of the Embayment. Accordingly, the crustal 
cross-section that was generalized from Figure VI. 2. B for the purposes 
of this study is given as the 4-body model shown in Figure VI. 2. C. The 
gravity analysis due to Cordell (1977) was used to project the charac- 
teristics of this generalized crustal cross-section north and south of 
the Y-S profile. 

Cordell (1977) corrected the smoothed positive Bouguer anomaly of 
the Embayment for the low-density sediments and observed the long con- 
tinuous positive anomaly with an amplitude of 15 to 45 mgal increasing 
southward that is illustrated in Figure VI. 2. D. The axis of this anomaly 
follows closely the Mississippi River northward beyond its confluence 
with the Ohio River into southern Illinois. The anomaly exhibits fairly 
uniform behavior south of the Y-P profile until about 33°N where it 
increases sharply, thus, suggesting that the crustal cross-section may 
be uniformly projected southward along the Mississippi River to approxi- 
mately 33°N. Northward decreasing gravity anomaly values in conjunction 
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with the northward tapering surface configuration of the Embayment 
suggest a commensurate northward tapering projection of the crustal 
cross-section along the Mississippi River into southern Illinois. 

Hence, to obtain the first-order, three-dimensional generalisation 
of the crustal structure of the Embayment used in this investigation, 
the crustal cross-section of Figure VI. 2. C was projected uniformly 
south of the Y-S profile and tapered uniformly northward as indicated 
in Figure VI. 1. The northern ends of the 4 bodies of this generalised 
model as projected onto the cross-section along the Y-S profile are 
given by the shaded regions of Figure VI. 2. C. 

To compute the potential field anomalies at 450 km elevation, 
each of the 4 bodies of this generalized model was represented by a 
Gauss-Legendre quadrature formula consisting of 128 equivalent point 
sources. The latitude and longitude limits of each body were repre- 
sented each by 8 point sources and the radial limits by 2 point sources. 
Pertinent body volume limits were interpolated from a set of body 
points that sampled the coordinates of the surface envelope for each 
body. The quadrature formulae were next evaluated and summed over a 
(21,13) observation grid spanning the region (260-280)°E, (33-45)°N 
and compared to observed gravity and magnetic anomaly data at 450 km 
elevation. 

D. Results and Discussion 

Free-air gravity anomaly values upward continued from the surface 
of the earth to an elevation (Z) of 450 km by equivalent point source 
inversion are illustrated for the study area in the stereographic 
equal-area polar (SEAP) projection in Figure VI. 3. A. The amplitude 
range (AR) of the data set is between 1.9 and -7.4 mgal and the amplitude 
mean (AM) is -3.6 mgal. These data exhibit a pronounced relative 
positive anomaly with slightly greater than 3 mgal of relative 
amplitude in the region of the Embayment. The gravity effect (Figure 
VI. 3. B) of the generalized 4-body model described above is roughly a 
4 mgal anomaly that, in general, corresponds to the observed 
data. The general degree of this correspondence for profiles along 35°N 
and 37°N is illustrated in Figures VI. 4 and 5, respectively. Here, 
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Figure VI. 3 Mississippi Embayment satellite-level comparisons between Gauss-Legendre quadrature modeled 
gravity and magnetic anomalies and upward continued l°-averaged free-air gravity and reduced 
to the pole POGO satellite magnetic anomaly data. Each anomaly field is plotted on a stereo 
graphic equal-area polar projection. 
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F i gu re V I , 4 Mississippi Embayment satellite— level profile comparisons between Gauss— Legendre quadrature 
modeled gravity anomaly and observed upward continued l°-averaged free-air gravity anomaly 
data along 35 N latitude. The shaded region indicates the location of the generalized 
Mississippi Embayment crustal model along the profile. 
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the profile peaks of the modeled gravity anomaly have been adjusted 
to the respective peak values of the observed profiles over the 
Embayment model to facilitate the comparison. 

A study of the profile comparisons shows increasing disparity 
between the modeled and observed data along the anomaly flanks away 
from the model. This may reflect contributions to the observed data 
from sources outside the study area that have not been modeled. None- 
theless, the general agreement between the modeled and observed data 
over the Embayment suggests that the observed gravity anomaly can be 
accounted for reasonably well at 450 km elevation by the 4-body model 
generalized in the previous section. 

POGO satellite magnetometer observations reduced to radial polari- 
zation using a normalization amplitude of 60,000 gamma by equivalent 
point source inversion are given for the study area in Figure VZ.3.C. 
These data show a prominent east-west magnetic high that is breached 
in the vicinity of the Embayment by a magnetic low. Inverse correla- 
tions of positive gravity and negative magnetic anomalies that char- 
acterize the Embayment also have been observed at satellite elevations 
for features such as the Yellowstone geothermal area where it is proposed 
that the net negative magnetization may be due to upward deflection of 
the Curie isotherm (von Frese et al . , 1979). In this regard, Mitchell 
et al., (1977) using travel time residuals of teleseismic P-waves 
also have speculated that the New Madrid seismic zone may be underlain 
by a mantle hot spot such as suggested for the Yellowstone region. 

To better resolve the characteristics of the magnetic anomaly for 
the Embayment region, the radially polarized data were high-pass 
filtered for anomaly wavelengths smaller than about 10°. The high- 
pass filtered data are illustrated in Figure VI. 3. D and show a negative 
anomaly of roughly -3 gamma over the Embayment. Wasilewski et al . , (1979) 
found that all of the information from analysis of medium to long- 
wavelength magnetic anomalies indicates that the sources probably are 
contained in the lower crust which, in general, may be substantially 
more magnetic than the upper crust. Typical estimates of deep crustal 
magnetization are on the order of 5x10 ^ emu/cm" (e.g.. Hall, 1974; 
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Shuey at al_. , 1973) . Wasilewski at al. * (1979) proposa that conditions 

for coharant regional magnatization are enhanced as crustal depth 

increases. Reman anca and thermal overprints are diminished* and 

viscous magnetisation and initial susceptibility are enhanced with 

o o 

increasing temperature especially within 100 C-150 C of the Curie point. 

The thickness of the crust within this thermal regime of the Curie 
point may be 5 to 20 km depending on the steepness of the geothermal 
gradient. Accordingly* they suggest that deep crustal magnetic sources 
probably are related to lateral variations of petrologic factors or 
Curie isotherm topography. 

Consideration of the foregoing remarks indicates that the most 
obvious deep crustal source for the observed magnetic anomaly is body 
#2 which also represents the major gravity source of the Embayment 
model. Austin and Keller (1979) propose that the combination of 
bodies #1 and #2 was formed as a manifestation of a mantle upwarp 
beneath the Embayment comprising of a mixture of crust and upper 
mantle material that cooled subsequently to form a block of high 
density material. Magnetic hypotheses which are consistent with this 
view include body #2 as a zone of negative magnetization contrast with 
respect to the lower crust due to depletion of magnetic minerals. 

Negative magnetization for body #2 also can result from temperatures 
which exceed the Curie point* although present heat flow data (Sass 
et al. , 1976) do not appear to warrant this hypothesis for the Embayment. 

Body #3 may represent an additional magnetic source for the Embay- 
ment assuming crustal magnetization increases with depth. However* 
the positive magnetic contribution of body #3 will be relatively weak 

i 

if the Curie isotherm depth i« about 40 km or more. Arguments for 

including bodies #1 and #4 in a magnetic model of the Embayment appear | 

i 

to be lacking, so that body #2 probably represents the primary source j 

r 

for the observed magnetic anomaly data if the Curie isotherm is at 

about 40 km of depth in the region of the Embayment. } 

i 

Accordingly, the magnetic anomaly due to body #2 was calculated in 

-3 3 ! ; 

Figure VI. 3. E using a magnetization contrast of -2.4x10 emu/cm . These 

i j 

results show that the anomaly amplitude observed for the region of the 1 \ 
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observed magnetic anomaly data along 35^ latituce. The shaded region m 
of the generalized Mississippi Embayment crustal model along the profile. 
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the Embayment at 450 km elevation can be well matched by a source ouch 
as body #2 located near the base of the crust with magnetic properties 
which correspond with the magnetization characteristics anticipated for 
the lower crust. Hence/ body #2 as determined for gravity modeling 
considerations represents at least a preliminary magnetic model for 
the Embayment . 

Substantial differences are apparent, however, between the spatial 
characteristics of the observed and modeled magnetic anomalies. Pro- 
file comparisons between Figures VI. 3. D and 3.E along 35°N and 37°N 
are given in Figures VI. 5 and 7, respectively, which illustrate 
this disparity. Here again, the modeled profile peaks are adjusted to 
the observed anomaly peaks over the Embayment model to facilitate com- 
parison. In general, these comparisons indicate that body #2 is mag- 
netically more extensive to the north and more restricted in longitude 
than it is gravitationally, although care must be taken with this inter- 
pretation because the observed data obviously contain the magnetic 
effects from sources outside the Embayment. Further refinements of 
the magnetic model are necessary and will be particularly warranted 
when the data are available from the current Magsat program (Langel, 

1979) to verify and further upgrade the POGO satellite magnetic anomalies 
for lithospheric applications. 

E. Conclusions 

A preliminary three-dimensional model of crustal structure for 
the Mississippi Embayment was generalized from analyses of surface 
wave dispersion (Austin and Keller, 1979) , seismic refraction data 
(McCamy and Meyer, 1966) and gravity data (Ervin and McGinnis, 1975; 
Cordell, 1977). The agreement between the gravity anomaly calculated 
for the 4-body model and the gravity data observed at 450 km elevation 
suggests that this model represents a reasonably valid generalization 
of the crustal density structure of the Embayment. The magnetic anomaly 
calculations show that the magnetic data observed at 450 km elevation 
can be modeled by a non-magnetic portion of the lower crust located 
along the axis of the Embayment that also corresponds to the major 
gravity source of the region. These results are consistent with the 
failed-rift hypothesis for the origin of the Mississippi Embayment. 
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This investigation further indicates that an aulacogen may be 
characterized at satellite elevations by observable positive gra.ity 
and negative magnetic anomalies. The primary source of both anomalies 
is the portion of the rift that defines a lateral variation in the 
physical properties of the lower crust. Finally, these results 
demonstrate the basic difficulty of separating pertinent anomalies 
at satellite elevations which perhaps is a major limitation to 
interpretation of satellite observed gravity and magnetic anomaly data 
However, the application of the modeling capacity demonstrated in 
this investigation will help to alleviate this problem. 
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VII. SUMMARY AND CONCLUSIONS 


Artificial, earth-orbiting satellites are making increasingly 
available consistent, regional-scale gravity and magnetic anomaly data 
which are largely free from the effects of near-surface geologic sources 
that otherwise tend to distort the signatures of deep, broad litho- 
spheric variations. Accordingly, satellite potential field obser- 
vations facilitate a unique characterization of the attributes of 
geologically significant regions measured in hundreds or even thou- 
sands of kilometers. These zones identified and characterized on a 
global basis can provide useful information for deciphering earth 
history including paleo and contemporary geodynamics, delineation of 
segments of the lithosphere into resource provinces, and for numerical 
modeling of lithospheric processes. 

To realize this promise, however, procedures are required to 
accurately and efficiently analyze potential field anomalies for 
lithospheric information directly in spherical coordinates. As demon- 
strated by this investigation, an especially powerful and efficient 
approach to generalized spherical earth processing of satellite potential 
field observations is to relate the gravity or magnetic anomalies to 
spherical distributions of equivalent point sources by methods of least 
squares matrix inversion. Linear transformations of the resultant 
equivalent source field permit generalized processing of the potential 
field anomalies in the spherical domain that is analogous to conven- 
tional geophysical processing methods routinely employed in the 
Cartesian plane. 

For example, spherical earth processing of gravity anomalies by 
equivalent point source inversion can yield corresponding geoidal 
deflections, vector anomaly components, anomaly continuations and 
differentiations at any point in source-free space. Similar elements 
for magnetic anomalies can be obtained in the spherical domain by 
equivalent point source inversion as well. In addition, equivalent 


156 


point source inversion can be used to differentially reduce regional 
magnetic anomalies to the radial pole. In general* the precision of 
these generalized spherical earth processing procedures is equivalent 
to the accuracy of the least squares fit achieved by the matrix inver- 
sion between the equivalent source field and the observed anomalies. 

Equivalent point source inversion was used to investigate regional 
correlations of surface free-air gravity and POGO satellite magnetic 
anomalies and regional heat flow and tectonic data for North America 
and adjacent marine areas. The combined analysis of these different 
classes of geophysical data is useful for improving the quality and 
uniqueness of their geologic interpretation, particularly in consid- 
eration of the general absence of regional ground-truth information. 

The qualitative analysis by visual spatial correlation indicates 
that regional gravity anomalies reflect the prevailing structural trends 
of both oceanic and continental terranes. When upward continued to 
450 km elevation, the gravity anomalies show strikingly direct cor- 
respondences to the regional heat flow pattern, whereas magnetic 
anomalies demonstrate limited correlation. 

Over continental terrane gravity and magnetic anomalies exhibit 
basically an inverse correlation at satellite elevations. One type 
of inverse correspondence is exemplified by the transcontinental mag- 
netic maxima and gravity minima that span the midcontinental region 
between the Anadarko Basin and Cincinnati Arch. Seismic evidence shows 
this region is characterized by enhanced crustal thickness, thus, 
suggesting that the regional combination of negative gravity and positive 
magnetic anomalies may reflect displacement of denser and less magnetic 
mantle material by an anomalously thick crust. 

The second type of inverse correlation is typified by the magnetic 
minimum and gravity maximum observed for the Yellowstone geothermal 
region. Here the regional inverse correspondence of positive gravity 
and negative magnetic anomalies may be related to increased heat flow 
triiich is likely to produce a region of relative net negative magnetization 
by deflecting the Curie isotherm to a shallower depth, as well as a 
thinned crust and inflated elevation that may not be isostatically com- 
pensated. 
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A more quantitative investigation of these regional correlations 
was developed utilizing Poisson's theorem in a moving-window linear 
regression analysis performed between suitable derivatives of the 
anomalous gravity and magnetic potentials at satellite elevations. 

These results are useful for establishing patterns of anomaly cor- 
relation and provide initial physical property constraints to facilitate 
further investigation of correlating anomalies. 

Over continental terrane, for example, the previously described 

inverse correlation of regional gravity and magnetic anomalies is well 

delineated by the regional patterns of the statistical correlation 

coefficients derived from the moving-window Poisson's analysis. In 

addition, these results indicate surprisingly large magnetizations 
-2 3 

(3 1 . 8 x 10 emu/cm ) for the transcontinental magnetic sources. 

Over the oceans, regional correlations of gravity and magnetic 
anomalies appear limited and poorly developed. The results of the moving- 
window Poisson's analysis suggest this observation may reflect decreased 
magnetic polarization of the oceanic crust relative to continental 
crust due to the influence of temperature and time on the remanent 
moments of crustal extrusive rocks. Statistically, however, oceanic 
regions are characterized primarily by direct regional anomaly correla- 
tions. This correspondence perhaps is best exemplified by the higher 
order gravity and magnetic anomaly derivatives of the Caribbean, eastern 
Gulf of Mexico and Atlantic Ocean. 

In general, the broad conclusions suggested by the correlation 
analyses developed during the course of this investigation indicate 
that long-wavelength magnetic anomalies as mapped by the POGO and more 
recently the Magsat programs provide a valuable complementary data set 
to gravity and other classes of geophysical data for investigating the 
megatectonic features of continents and oceans. 

The capacity to directly model potential field anomalies observed 
in spherical coordinates was investigated as an essential element to 
complement the spherical earth processing and correlation procedures 
described above. Accordingly, a technique was developed and demon- 
strated which employs Gauss-Legendre quadrature integration to numer- 
ically model gravity and magnetic anomaly potentials and their derivatives 
in the spherical domain due to arbitrarily shaped geological bodies 
■< with arbitrary physical properties. 
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The modeling procedure involves representing the anomalous body 
by a distribution of equivalent point sources which is determined 
directly from the volume limits of the body. Calculating and summing 
at each observation point the effects of the equivalent point sources 
located at these selected coordinates within the body's volume gives 
an accurate estimate for the potential field anomaly of the arbitrary 
body. The method is especially well suited for satellite gravity and 
magnetic anomaly modeling because the efficiency and accuracy of the 
application increases with increasing distance between source and 
observation points. 

The Gauss-Legendre quadrature modeling procedure was used to 
investigate a 4-body model of northern Mississippi Embayment crustal 
structure that was generalized from analyses of surface wave disper- 
sion, seismic refraction and surface Bouguer gravity anomaly data. 

The modeled gravity anomaly signature compares favorably to the 
relative positive free-air gravity anomalies upward continued to 450 km 
elevation, thus, suggesting that the 4-body model is a geologically 
reasonable representation of the density structure of the Embayment. 

Magnetic modeling shows that the pronounce minimum observed in 
the POGO satellite magnetic anomaly data differentially reduced to 
radial polarization can be accounted for by a block of non-magnetic 
material located at the base of the crust along the axis of the Embay- 
ment. This particular block of lower crust also corresponds to the 
primary gravity source for the region according to the gravity modeling 
considerations . 

The results of the modeling application are consistent with the 
failed-rift (aulacogen) hypothesis for the origin of the Mississippi 
Embayment. More generally, these results suggest that an aulacogen 
may be characterized at satellite elevations by observable positive 
gravity and negative magnetic anomalies, due principally to the rift 
component which defines a non-magnetic block of high-density material 
within the lower crust. 

In contrast to the widespread use of potential field anomalies 
referenced in the Cartesian plane for geologic analysis, there has 
been in the past a general inability to routinely analyze gravity 
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and magnetic anomalies in spherical coordinates for subsurface infor- 
mation that has limited the utility of satellite potential field obser- 
vations in lithospheric studies. During the course of this investiga- 
tion, however, procedures have been developed and verified for generalized 
processing, correlation analysis and modeling of potential field 
anomalies in the spherical domain that are analogous to the essential 
geophysical methods of interpretation commonly employed in conventional 
flat earth applications. 

In general, then, these spherical earth techniques can be used to 
routinely perform on spherically registered potential field anomalies 
nearly any type of geophysical analysis that finds widespread use in 
the flat earth domain. Accordingly, it is concluded that the ultimate 
utility of satellite potential field observations for lithospheric 
analysis now is basically a technological question of how accurately 
lithospheric potential field signals can be measured and reduced by 
satellite programs. Clearly, the results and procedures of this inves- 
tigation can play a significant role in helping to resolve this survey 
design problem. 

In summation, the results and procedures developed during the 
course of this investigation have been shown to represent a powerful 
and efficient approach to the geological analysis of gravity and magnetic 
anomalies in the spherical domain. Accordingly, it is concluded that 
the methods and results of this study have widespread application in 
the analysis and design of regional-scale gravity and magnetic surveys 
for lithospheric investigation. 
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